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This  work  describes  the  development  of  a  computer  based 
model  that  would  allow  for  determination  of  the  transient 
response  characteristics  of  jet  vanes  of  any  size.  The  model 
used  a  thermal  lump  approach  method,  considering  only 
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data  obtained  from  test  firings  of  a  rocket  motor.  The  full 
scale  modeling  results  were  compared  to  previous  quarter  scale 
models  in  am  atten^t  to  determine  the  applicability  of  the 
quarter  scale  results  to  full  scale  vanes.  It  was  determined 
that  the  quarter  scale  model  did  not  provide  an  accurate 
representation  of  the  heat  transfer  process  in  larger  scale 
vanes,  although  the  full  scale  model  provided  a  sufficient 
representation  of  the  thermal  response  of  the  jet  vane. 
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I.  INTRODUCTION 


The  work  reported  in  this  document  is  further  research  in 
the  study  of  Thrust  Vector  Control  (TVC)  system  used  in 
missile  trajectory  guidance  during  initial  phases  of  a 
missiles  flight  [Ref.  1].  TVC  utilizes  a  system  which 
deflects  exhaust  thrust  to  control  the  missile  during  launch 
and  initial  flight  stages.  For  this  phase  of  flight  the  vane 
is  actuated  from  a  stowed  position  and  inserted  into  the 
exhaust  of  the  rocket  nozzle  to  provide  the  required  control. 
Figure  1.1  illustrates  the  stowed  and  active  positions  of  the 
TVC  vane. 


Figure  1.1  STARS  TVC  Sto«#age  Concept 


1 


The  vane  operates  in  an  extremely  harsh  thermal 
environment.  In  order  to  use  adequate  materials  for  vane 
design  the  heat  transfer  characteristics  of  a  selected 
material  operating  in  such  an  environment  needs  to  be 
resolved.  In  order  to  determine  the  best  performance 
characteristics  of  different  materials  a  mathematical  model 
needs  to  be  constructed.  The  model  needs  to  contain  both  the 
known  physical  properties  of  the  material  and  the  unknown 
properties  to  be  determined.  It  will  utilize  the  experimental 
temperature  time  histories  from  rocket  firings  and  adjust  the 
unkno%ms  physical  properties  until  the  predicted  temperature 
time  history  and  the  experimental  data  are  close  in  a  least 
squares  sense. 

Past  work  in  the  system  identification  process  has  been 
conducted  using  a  quarter  scale  model  to  determine  the 
convective  heat  transfer  coefficients  and  to  verify  the 
applicability  of  the  quarter  scale  results  in  predicting  the 
heat  transfer  characteristics  of  a  full  scale  protoype.  The 
idea  behind  this  report  is  to  enlarge  the  model  to  full  scale 
in  order  to  obtain  the  heat  transfer  properties  and  to  compare 
these  results  with  those  obtained  from  the  quarter  scale 
modeling  process. 
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II.  BACKGROUND 


A.  LUMPED  PARAMETER  MODELING 
1.  Vane  Mbdeling  Concept 

The  jet  vane  that  is  part  of  the  TVC  system  concept 
will  be  used  as  a  deflection  device  to  direct  the  thrust  of 
the  rocket  exhaust.  The  vane  will  be  encountering  effects 
from  such  things  as  con^lex  thermal  changes,  supersonic 
particle  laden  stream  flow,  shock  waves  impinging  on  the  vane 
and  phase  variations  resulting  from  vane  ablation  [Ref.  2:p. 
5] .  The  idea  of  vane  modeling  is  a  process  by  which  some  or 
all  of  these  processes  may  be  modeled  in  order  to  find  the 
material  that  provides  the  durability  and  the  necessary 
performance  characteristics.  Due  to  the  mathematical 
complexity  and  the  wide  number  of  variables  that  may  effect 
vane  performance  it  is  necessary  to  limit  the  model  to  a 
smaller  scale.  The  main  emphasis  of  this  report  is  to  limit 
the  model  to  determining  the  thermal  convective 
characteristics.  This  modeling  concept  involves  development 
of  a  thermal  model  taking  into  account  the  conductive  and 
convective  effects  of  the  material  selected.  Radiation 
effects  will  not  be  considered  in  the  development  of  this 
model.  The  model  uses  the  concept  that  the  transient  thermal 
response  of  the  vane  can  be  developed  into  a  set  of  equations 
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that  describe  the  heat  transfer  process  occurring  over  a 
distinct  period  of  time. 

2.  Lumped  Parasoeter  Approach 

The  thermal  modeling  approach  taken  will  be  to  break 
the  vane  down  into  distinct  thermal  lunps  (Ref.  2;p.  101. 
This  approach  lumps  the  heat  transfer  parameters  into  distinct 
nodes .  The  geometric  data  and  physical  properties  of  the  vane 
material  at  are  then  used  to  determine  the  conductive 
resistance  and  capacitances  at  each  node.  The  convective 
resistances  between  the  exhaust  plume  and  the  material  are 
essentially  unknown  and  will  be  determined  using  the  modeling 
process.  The  lumped  parameter  approach  is  essentially  a 
division  of  the  vane  into  critical  locations <  usually  chosen 
to  correspond  to  a  thermocouple  attachment  point  on  the 
experimental  model <  in  order  to  simply  the  modeling  process 
(Ref.  2:p.  12].  This  simplification  enables  the  formulation 
of  energy  balance  equations*  mentioned  previously*  that 
describe  the  transient  thermal  response  of  the  individual 
nodes  in  the  vane  (Ref.  l:p.  4].  The  energy  balance  equations 
take  into  account  thermal  energy  .ito  and  out  of  the  vane 
nodes  and  the  temperature  of  the  vane  accounts  for  the  energy 
stored  at  the  nodes. 

3.  Parametric  System  Identification 

The  thermal  model  will  be  computer-based  FORTRAN  code 
using  Parametric  System  Identification  (PSD  (Ref.  3:p.  2]. 
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PSI  is  a  procedure  in  which  the  model  parameters  generate  a 
temperature  time  history  that  predicts  the  experimental 
tenperature  time  history  [Ref.  2:p.  6]  .  The  process  involves 
adjustment  of  the  unknown  parameters  until  a  desired  agreement 
is  reached  between  experimental  and  model  temperature 
profiles.  The  PSI  process  involves  adjustment  of  the  unknown 
values  for  the  convective  coefficients  until  the  model 
temperature  time  histories  of  selected  nodes  are  matched  with 
the  experimental  temperature  time  histories  of  the 
corresponding  nodes  on  the  actual  vane.  The  PSI  process 
brings  together  the  mathematical  model  and  the  experimental 
data.  The  accuracy  of  the  match  between  model  and 
experimental  data  is  limited  to  the  complexity  of  the  energy 
equations  used  to  describe  the  model  itself.  The  process  can 
be  a  self  adjusting  one  in  that  the  model  can  be  used  to 
design  the  vane  for  experimental  testing  and  then  the 
resulting  experimental  data  can  be  used  to  update  the  original 
model  (Ref.  2:p.  6]. 

B.  B^slc  mosLim  processes 

1.  Ecruatioo  Porxnulation 

The  modeling  of  the  TVC  vane  involves  generating  the 
energy  balance  equations  for  the  lumped  parameter  model.  The 
energy  equations  are  formulated  by  considering  the  flow  of 
energy  into  and  out  of  a  node.  Figure  2.1  illustrates  the 
energy  flow  concept  of  a  vane  node. 
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Figure  2.1  Thezinal  Node  Energy  Flow 
The  energy  balance  formulation  for  Figure  2.1  can  be 
constructed  by  considering  the  flow  of  energy  into  and  out  of 
the  node.  This  energy  balance  is  represented  by, 


f  .JL-JL-JL.JL 

'  CiRi  C^Ro  ^iRo 


(1) 


where  the  value  of  Tj>Ti>To.  Starting  with  the  determination 
of  the  resistances  and  capacitances  at  each  individual  node 
from  the  geometric  data  and  the  physical  properties  of  the 
material  is  required.  Prior  to  understanding  the  effects  of 
sub-scale  model  the  conductive  and  convective  parameters  first 
need  to  be  addressed.  The  conductive  resistance  is  found  by, 

where  L  is  the  conductive  path  length,  \  is  the  cross 
sectional  area  of  the  conductive  path  and  K  is  the  thermal 
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conductivity  [Ref.  3:p.  7].  Likewise,  the  convective 

resistance  can  be  found  ly. 


{3) 

where  A, 

is  the  convective  surface 

area 

and  h  is 

the 

convective 

heat  transfer  coefficient 

[Ref. 

3:p.  7]. 

The 

thermal  capacitance  is  found  by, 

(4) 

where  q  is  the  material  density,  Cp  is  the  specific  heat  and 
V  is  the  volume  [Ref.  3:p.  7].  It  can  be  seen  that  these 
values  are  dependent  on  not  only  the  material  in  which  the 
vans  is  constructed  but  also  the  actual  size  of  the  vane  being 
modelled.  All  geometric  data  is  given  in  full  scale  and  the 
physical  properties  are  assumed  constant  no  matter  what  size 
model  is  being  considered. 

The  vane  model  is  generated  by  formulating  the  energy 
balance  equations  for  each  individual  node  of  the  vane.  The 
equations  will  be  dependent  on  the  number  of  nodes  used  to 
model  the  heat  transfer  process  occurring  in  the  vane.  To 
develop  such  a  set  of  equations  the  nodes  locations  need  to  be 
established  on  the  vans.  The  nodes  chosen  should  correspond 
to  a  center  of  mass  or  the  major  sections  of  the  object  being 
modeled,  in  this  modeling  process  the  vane  is  broken  down 
into  four  distinct  sections.  Figure  2.2  illustrates  a  vane 
with  a  four  node  configuration  [Ref.  l;p.  6). 
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Figure  2.2  Four  Node  Vane  Model 
This  four  node  select: ion  is  by  no  means  a  standard  which  must 
be  adhered  to.  Multiple  nodes  may  be  selected  depending  on 
the  complexity  that  is  desired. 

Xn  Figure  2.2  the  vane  is  broken  down  into  four  distinct 
nodes.  Node  one  is  for  modeling  the  tip,  node  two  for  the 
fin,  node  three  for  the  shaft  and  node  four  for  modeling  the 
mount  (Ref.  l:p.  5).  The  energy  balance  equations  are 
generated  by  considering  the  heat  transfer  into  and  out  of 
each  individual  node.  Nodes  one  and  two,  two  and  three,  three 
and  four,  are  connected  by  conductive  paths  with  a  conductive 
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resistance  linking  the  nodes  (Ref.  3:p.  4].  The  designation 
Rij  is  used  to  denote  the  conductive  resistance  between  nodes 
i  and  j .  Nodes  one  and  two  are  linked  convectively  to  two 
driving  temperatures  TRl  and  TR2 .  These  are  known  as  the 
recovury  temperatures  and  are  the  source  of  the  heating 
process  [Ref.  3:p.  4].  TRl  corresponds  to  the  stagnation  flow 
while  'rR2  represents  the  tenperature  of  the  free  stream  flow. 
Tue  convective  resistances  are  represent  by  the  designation 
Rpi,  where  the  i  represents  the  node  that  the  convective 
temperature  is  driving.  Each  node  contains  a  thermal 
capacitance  represented  by  Ci  where  i  is  the  node 
representation.  The  generation  of  the  energy  balance 
equations  involves  looking  at  each  individual  node  and 
applying  the  law  of  conservation  of  energy  [Ref.  3:p.  4). 
The  conservation  of  energy  for  each  node  involves  energy 
transfer  in  and  out  of  the  node.  That  is  the  rate  of  change 
of  energy  at  a  node  is  equal  to  the  er.arny  into  the  node  minus 
the  energy  out  of  the  node.  The  remaining  energy  at  the  node 
is  what  causes  the  temperature  change  at  that  particular 
location.  The  node  is  a  point  loca^-ion  and  the  transient 
thermal  heating  is  being  considered. 

/^splication  of  this  idea  leads  to  the  following  equations 
[Ref .  3 :p.  4) ; 


^i^Pi  ^1^2 


(5) 
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(6) 


In  order  to  make  conputer  coding  easier  a  sin^lification  of 
the  above  equations  is  desired.  The  following  characteristic 
rate  coefficients  are  defined  (Ref.  3:p.  9]: 


1 

(9) 

«3a“ 

1 

€^^23 

’*  CA. 

(10) 

^40“ 

1 

jb,,  . 

(11) 

A  further  simplification  can  be  made  by  combining  coefficients 
for  the  same  temperature  at  a  particular  node.  This 
simplification  allows  each  node  temperature  T|  to  have  a 
single  coefficient  making  it  possible  to  put  the  energy 
balance  equations  into  matrix  form.  This  is  done  as  follows 
(Ref.  3;p.  5];  The  energy  balance  equations  now  become  (Ref. 
3  s  p  •  5  ] 
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^33  “^32  ■*■^34  (13) 

(16) 

tl=a43  ^3-344^4  (17) 


These  equations  may  be  put  in  matrix  form  to  give  a  clearer 
understanding  of  there  use  in  the  computer  coding 
(Ref .  3 :p,  5] , 

■  r,  ]  -an  On  0  ^  1  f  1  \  0  0  0  ]  T  Tm  ' 

Tj  _  an  -On  aja  0  ^2  ,  0  622  0  0  Tm 

T3  0  a3j  —033  O34  Ta  0  0  0  0  0 

ti,  L  0  0  043  -044  J  I  T4  .  .  0  0  0  0 .  .  0  . 

This  is  a  compact  notation  commonly  known  as  a  state  space 
equation  and  may  be  written  as, 

jic°AJC'^Bu  (19) 

Equations  of  this  form  are  common  and  a  wide  variety  of 
package  are  available  for  solving  linear  differential 
equations . 
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2 .  Previous  Slodels 


Work  on  the  jet  vane  thermal  model  was  started  by  Nunn 
and  Kelleher  in  1986  at  the  Naval  Postgraduate  School  [Ref. 
5]  .  Further  research  into  the  model  was  conducted  by  Nunn 
[Ref.  2],  Hatzenbuehler  [Ref.  4]  ,  Reno  [Ref.  1]  and  Driels 
[Ref.  3].  Prior  work  has  concentrated  on  development  of  a 
quarter  scale  model  and  its  applicability  to  full  scale  vanes. 
Reno  used  both  3  and  4  node  approaches  in  modeling  the  heat 
transfer  process  of  the  vane.  The  experimental  data  used  for 
the  comparison  process,  on  the  quarter  scale  model,  was  from 
the  shaft  and  mount  locations  [Ref.  l:p.  7] .  The  points  were 
chosen  were  due  to  the  limited  size  of  the  experimental  veuie 
for  thermocouple  attachment.  Those  points  being  the  back  of 
the  shaft  and  mount.  The  4  node  vane  model  in  figure  2.1 
illustrates  the  set-up  used  by  Reno  along  with  the  energy 
balance  equations.  The  unknovm  variables  that  needed  to  be 
determined  were  Rp2#  ^34*  K4Q  and  C4  [Ref.  l:p.  10]  .  In  the 
simplified  energy  balance  equations  these  coincide  to  a34,  a43, 
a4Q  and  b22<  The  geometric  data  used  for  the  full  scale  vane 
is  given  in  Table  1  [Ref.  3:p.  8]. 

TABLE  1 


tip  to  vane  to  shaft 

V,  cm^ 

A^,  cm^ 

L,  cm 

2.6  52  23 

5.9  5.2 

5,0  6.0 
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The  physical  data  for  the  vane,  which  was  a  tungsten  matrix 
impregnated  with  copper,  is  p  =  18310  kg/m^  K  =  173  W/m-K, 
and  Cp  =  146  J/kg*K  [Ref.  3;p.  7].  The  recovery  tenperatures 
used  were  TRl  =  2360  K  and  TR2  =  2260  K.  The  input  thrust 
profile  used  was  simulated  as  a  ramp  function.  The  recovery 
temperatures,  TRl  and  TR2,  contained  in  the  input  vector  u,  of 
the  matrix  form,  drive  the  system,  and  are  the  product  of  the 
steady  state  values  and  the  thrust  function.  Figure  2.3 
illustrates  this  functional  input  [Ref.  3]. 


Figure  2.3  Rasp  Function  for  Recovery  Tenperatures 
Reno  utilized  a  personal  computer  with  a  program  System  Build 
and  System  Identification  developed  by  Integrated  Systems, 
Inc.  [Ref.  l:p.  10]  .  The  identification  process  used  MAXLIKE 
which  is  a  function  of  MATRIX  X  [Ref.  l:p.  10].  This  program 
allowed  the  four  unknown  parameters  to  vary  until  the 
temperature  profiles  generated  matched  those  obtained  from  the 
experimental  firings.  The  following  results  were  obtained: 
bjj  o  0.2718  s’‘,  aj4  =  0.299  s'‘,  a*,  =  0.1322  s‘‘  and 
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340  =  0.1018  s'^  [Ref.  I:p.l0].  These  parameters  correspond  to 
the  values  of  [Ref.  l:p.  11]:  R„  =  1.67  K/W,  R34  =  3.33  K/w, 
R40  =  4.33  K/W,  C4  =  2.27  J/K.  There  was  a  good  match  between 
the  generated  and  experimental  temperature  time  histories  of 
the  mount  and  the  shaft.  Figure  2.4  demonstrates  the  close 
proximity  between  the  shaft  and  mount  profiles  (Ref.  l:p.  12]  . 


TtNF.  (K) 

Figure  2.4  Quarter  Scale  Tenperature  Profiles 
Determination  of  the  unknown  paramet<  s  is  based  upon  a 
temperature  profile  match  between  the  mount  and  shaft.  The 
temperature  profiles  ten^eraturcs  and  Tj  are  the 

mathematical  equation  profiles  and  have  no  experimental  data 
to  be  matched  to  in  the  quarter  scale  modeling  process. 
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Work  done  by  both  Reno  (Ref.  1]  and  Driels  [Ref.  3] 
indicated  that  although  the  convective  heat  transfer 
coefficient  at  the  tip  node  was  an  unknown  variable,  in  the 
quarter  scale  modeling  process,  it  had  minimal  effect  on  the 
outcome  of  the  previously  identified  unknown  values  and  thus 
bii  was  held  constant.  Table  2  illustrates  the  variation  in 
the  four  unknowns  for  various  values  of  bn  [Ref.  3:p.  12]. 

TABLE  2 


:l9i 

a34 

^43 

^40 

^22  1 

5.35 

0.483 

0.152 

-0.167 

0.143  1 

10.0 

0.447 

0.152 

-0.168 

0.162 

20.0 

0.44 

0.152 

-0.169 

0.174 

1.0 

0.51 

0.153 

-0.164 

0.056 

0.5 

0.47 

0.153 

-0.163 

0.019 

The  results  reported  in  Table  2  are  those  obtained  by  Driels 
[Ref.  3]  and  do  not  exactly  match  the  results  obtained  by  Reno 
[Ref.  1].  The  tert®)erature  profiles  obtained  by  Driels 
coincide  with  those  displayed  in  Figure  2.4. 


C.  MODEL  SCALING 

1.  Purpose  of  Quarter  Scale  Model 

The  quarter  scale  model  was  developed  to  predict  the 
behavior  of  the  heat  transfer  pareimeters  of  a  full  scale 
prototype.  The  quarter  scale  model  quantitatively  provides 
for  testing  on  a  smaller  scale  which  in  the  long  term  is  more 
cost  effective  than  constructing  one  of  full  scale.  This 
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savings  come  from  the  reduced  size  of  the  rocket  motor,  the 
reduced  use  of  data  acquisition  equipment,  environmental 
requirements  that  have  to  be  met  when  firing  full  scale  that 
do  not  apply  to  the  quarter  scale,  manufacturing  of  quarter 
scale  vanes,  material  requirements,  etc.  The  quarter  scale 
modeling  process,  if  adequate,  could  provide  a  quick,  simple 
form  of  modeling  for  a  larger  vane  without  the  added 
difficulties . 

2.  Quarter  to  Full  Scale  Model  Predictions 

It  has  been  demonstrated  that  a  sub-scale  model  can  be 
generated  and  subsequent  experimental  data  obtained  in  order 
to  predict  the  heat  transfer  parameters  of  a  given  material . 
A  set  of  convective  heat  transfer  coefficients  can  be 
determined  from  the  results  of  such  models.  The  accuracy  of 
these  results  and  there  applicability  to  a  full  scale  vane  has 
yet  to  be  determined.  Direct  scaling  of  the  geometric 
properties  is  achieved  and  the  physical  properties  such  as 
density,  thermal  conductivity  and  specific  are  assumed 
constant  for  both  sub-scale  and  full  scale  vanes.  The 
difficulty  arises  in  the  effects  of  the  jet  plume  on  the 
different  sizes  of  vane.  As  mentioned  earlier  the  conductive 
properties  rely  on  both  physical  and  geometric  data. 
Variation  in  the  geometric  size  would  have  a  profound  impact 
on  the  transient  heat  transfer  process  occurring  not  only 
inside  the  vane  but  due  to  the  larger  surface  area  the  effects 
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of  the  convective  heat  transfer  process.  A  direct  scaling  up 
of  the  quarter  scale  model  to  full  scale  may  be  inadequate  to 
properly  model  the  effects  of  a  large  vane.  Increases  in  the 
rocket  chamber  temperature  thus  an  increased  recovery 
temperature  may  also  influence  the  energy  transfer  process. 
The  purpose  of  this  report  is  to  develop  a  full  scale  model, 
compare  the  results  to  those  obtained  by  the  quarter  scale  and 
determine  if  the  results  obtained  from  the  sub-scale  tests 
reflect  the  actual  heat  transfer  process  of  a  full  scale  vane. 

3 .  Scaling  Analysis 

To  scale  the  resistances  and  capacitances  a  reduction 
in  the  linear  dimensions  of  the  vane  must  be  calculated  in 
order  that  the  model  reflects  the  actual  size  of  the 
manufactured  vane  in  the  experimental  phase  [Ref.  3:p.  8] ,  A 
scale  factor  SF,  where  SF<1,  is  introduced  to  reduce  the 
geometric  dimensions  of  the  vane  and  subsequently  altering  the 
resistances  and  capacitances  [Ref.  3:p.  8].  The  scale  values 
can  be  determined  by  manipulating  equations  (2),  (3)  and  (4). 
When  manipulating  these  equations  the  physical  properties  such 
as  the  thermal  conductance  k,  the  specific  heat  Cp,  and  the 
density  q  remain  constant  and  do  not  change.  The  equation 
manipulation  results  in, 


'  JCA. 


L-SF  ^ 


_L _ 3^ 

KA,SF 


(20) 
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(21) 


where  Rf  is  the  full  scale  resistance  and  R,  is  the  sub-scale 
resistance.  An  important  point  to  note  is  that  the  scale 
factor  SF  being  less  than  one  causes  the  resistance  to 
increase  during  the  scaling  process.  An  exanple  of  this  is 
the  quarter  scale  models  previously  constructed  where  the 
scaled  resistance  R,  is  equal  to  the  full  scale  divided  by 
0.25  (1/4  scale  modeling).  This  results  in  the  scaled 
resistance  being  four  times  the  full  scale  value.  The  thermal 
capacitance  has  the  same  linear  dimension  scaling  as  the 
resistance  did,  which  is  given  by, 

Cf^pCpV  (22) 

C^*>Cf’SF^  (23) 

where  C(  is  the  full  scale  capacitance  and  as  indicated  is 
the  sub-scale  capacitance  (Ref.  3:p.  8).  There  is  no  direct 
scaling  of  the  convective  resistance  in  the  model  due  to  the 
fact  that  the  heat  transfer  coefficient  h  is  unknown.  The 
only  scaling  used  on  the  convective  values  is  if  convergence 
has  been  achieved  and  the  value  of  the  convective  resistance 
is  needed.  The  scaling  process  for  the  convective  resistance 
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is  similar  to  those  previously  stated  and  is  given  by  [Ref. 
3:p.  81, 


Although  there  is  an  increase  in  the  resistance  due  to  scaling 
and  a  decrease  in  the  capacitance,  the  CCTtbined  effect  of 
these  two  values  must  be  emphasized.  The  product  of  RC  is  a 
time  constant  for  a  particular  driving  temperature  [Ref.  l:p. 
7] .  Scaling  of  the  values  of  R  and  C  results  in  a  change  in 
the  time  constants  for  the  energy  balance  equations  from  the 
full  scale  to  the  sub-scale.  The  resultant  change  is 
demonstrated  in  the  following  formula; 

-  R^C,°RtCf'SF^  (25) 
SF 

The  time  constant  for  the  sub-scale  values  decreases  as  the 
product  of  the  full  scale  time  constant  times  the  square  of 
the  scale  factor.  It  will  be  seen  that  in  formulating  the 
energy  balance  equations  for  the  mathematical  model  that  the 
scaling  of  the  time  constants  will  have  a  profound  effect  on 
the  conductive  and  convective  values. 


19 


III.  FULL  SCALE  M}DELING 


A.  SIX  NODS  MODEL  APPROACH 

1 .  Concept 

As  discussed  previously,  prior  quarter  scale  modeling 
approaches  involved  using  three  and  four  node  models. 
Placement  of  the  node  points  in  the  model  were  dependent  on 
the  center  of  mass  at  determined  critical  locations  and  the 
availability  of  experimental  data  from  a  quarter  scale  test 
firing.  The  experimental  data  coming  from  the  shaft  and  mount 
locations.  These  temperature  profiles  were  matched  by  the 
matheinatical  model  and  the  results  used  in  determining  the 
convective  heat  transfer  coefficient  for  the  fin  section,  the 
unknown  resistance  between  the  mount  and  shaft,  and  the 
unknown  capacitance  of  the  mount .  A  single  value  of 
convective  heat  transfer  coefficient  was  obtained  from  the 
quarter  scale  tests  but  due  to  the  lack  of  experimental 
temperature  profiles  from  that  particular  location  the 
applicability  of  this  value  to  a  larger  vane  was  in  question. 

The  basis  for  increasing  the  nuitdaer  of  nodes  when 
attenpting  to  model  the  full  scale  vane  was  two  fold;  1)  more 
node  locations  would  give  a  more  accurate  representation  of 
the  actual  heat  transfer  process  going  on  in  the  vane  and  2) 
these  added  node  locations  had  experimental  data  available 
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which  allowed  for  a  more  diverse  set  of  temperature  time 
profiles  that  could  be  used  in  constructing  the  mathematical 
model.  The  21  locations  that  experimental  data  was  taken  from 
allows  a  greater  flexibility  in  the  selection  of  temperature 
profiles  to  be  matched.  Although  Nunn  [Ref.  2:p  12] 
recommended  that  for  simplicity  and  efficiency  the  minimum 
number  of  nodes  should  be  used,  to  accurately  model  the  heat 
transfer  process  of  a  full  scale  vane  a  larger  number  of  nodes 
was  needed.  The  modification  was  made  by  adding  two  more  node 
locations  to  the  four  node  model  in  the  fin  area.  Figure  3.1 
shows  the  six  node  model  design. 


TR2 

0 


Figure  3.1  Six  Kode  Vane  Model 
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2.  Geometric  Modification 

The  addition  of  two  nodes  to  the  fin  in  the  six  node 
model  require  a  modification  of  the  geometric  properties  in 
order  that  the  proper  volumes,  areas  and  lengths  are  used  in 
calculating  the  resistances  and  capacitances.  These  new 
values  are  given  in  Table  3, 

TABLE  3 


tip  to  node2  to  node3  to  node4  to  shaft 

V,  cm^ 

A,,  cm^ 

L,  cm 

2.6  12.86  17.15  21.44  23.0 

5.6  5.9  6.2  5.2 

2.774  2.774  2.774  5.0 

The  values  are  based  on  those  obtained  from  Table  1  [Ref.  3:p. 
8]  and  the  four  node  model  design.  The  parameters  for  the 
node  volumes  were  obtained  by  taking  the  vane  volume  and 
estimating  a  20%  change  in  the  volume  from  the  rear  of  the 
vane,  node  four  section,  to  the  forward  part  of  the  vane,  node 
two  section,  with  the  rear  being  the  largest  and  the  basis  for 
the  calculations.  The  cross  sectional  areas  were  determined 
in  the  same  manner  with  the  area  at  node  two  being  considered 
to  be  20%  less  than  node  three  and  node  four  cross  sectional 
area  to  be  20%  greater  than  that  of  node  three.  The  cross 
sectional  area  at  node  three  was  known  because  it  corresponds 
to  node  two  of  the  quarter  scale  four  node  model .  The  quarter 
scale  model  parameters  were  given  in  full  scale  terms  and 
scaled  down  in  the  program  formulation. 
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3.  Energy  Balance  Equation  Foirmulation 

Formulation  of  the  energy  balance  equations  is  done  by 
applying  the  law  of  conservation  of  energy  to  the  six  nodes  of 
the  vane.  Application  of  this  law  leads  to  the  following, 
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^5^45 

^S'^se 

^5^56 

(26) 


(27) 


(28) 


(29) 


(30) 


‘6  . 


^6^6 


(3X) 


As  discussed  in  Chapter  II  this  series  of  equations  can  be 
simplified  by  defining  a  set  of  characteristic  rate 
coefficients.  These  coefficients  relate  the  properties  of  the 
vane  to  the  time  constants.  The  simplification  involves  both 
a  conductive  coefficient  labeled  a^^j  and  a  convective 
coefficient  Icibeled  b^j.  This  simplification  leads  to, 
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These  simplifications  now  allow  the  energy  balance  equations 
to  be  put  into  matrix  form.  The  modification  provides  an 
easier  approach  to  the  program  formulation.  The  matrix 
equation  description  was  described  in  Chapter  II  and  with  the 
application  of  the  characteristic  rate  coefficients  the 
lumping  of  the  common  ten^erature  coefficients  in  each  energy 
balance  formulation  the  following  equations  can  be  used  to 
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generate  the  FORTRAN  code  used  to  solve  these  differential 
equations.  The  equations  become. 


(39) 

^2  ~^22 ^2  ■*■^23  ^3  '*'^2 ^S2 

(40) 

^3  “^32  ^2  “®33  ^3  '*‘^34  ^4  "'■■^3 

(41) 

^4  "^43  ^3  "^44  ^4  ■‘■^45  ^5  '*■■^4  ^R2 

(42) 

%  ~^54  ^4  ”®55  ^5  ■*‘®56  ^6 

(43) 

(44) 

4 .  Program  Formulation 

The  energy  balance  equations  are  a  set  of  liner 
ordinary  differential  equations  with  some  of  the  rate 
coefficients  as  unknowns.  The  process  to  evaluate  these 
unknowns  is  to  start  with  a  set  of  initial  conditions  for 
these  unknowns,  solve  the  differential  equations  over  a  period 
of  time  to  obtain  a  tenperature  time  history  for  each  node. 
We  then  con^are  the  experimental  data  from  a  set  of 
corresponding  nodes  and  the  results  obtain  from  the  model 
progreun  and  calculate  the  difference  of  the  sum  of  the  squares 
between  the  two  temperature  profiles.  The  programs  optimizer 
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adjusts  the  unknovm  variables  and  solves  the  differential 
equations  again  to  obtain  a  new  temperature  time  history. 
This  process  is  illustrated  in  Figure  3.2 


Figure  3.2  Block  Diagram  of  Modeling  Program 
The  blocks  labeled  Pre.dat,  Trans. for,  TV. mat, TM. mat  and 
Pro. m  will  be  discussed  in  the  following  section.  The  program 
PIO  reads  the  experimental  data  from  the  Datam.dat  file.  It 
then  reads  in  the  full  scale  geometric  properties  of  the  vane 
and  the  physical  properties  of  the  material  the  vane  is 
constructed  of  from  Xnput.dat.  It  calculates  the  resistance 
and  capacitances  for  the  nodes  and  the  rate  coefficients  a^j. 
The  initial  conditions  are  set  for  the  unknown  parameters  and 
PID  then  calls  the  optimizer  subroutine.  The  optimizer 
subroutine  intern  calls  another  subroutine  called  TEMP  which 
calls  the  needed  equation  solver  and  function  program  that 
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provides  the  differential  equations.  DIVPRK  (Double-precision 
Initial  Value  Problem  Runge-Kutta)  is  an  IMSL  routine  used  to 
solve  linear  ordinary  differential  equations  with  the  given 
set  of  initial  conditions.  The  subroutine  that  supplies  the 
differential  equations  is  called  FCN.  The  values  of  the 
driving  temperatures  are  supplied  by  the  main  program  PID  and 
are  common  to  all  the  subroutines.  The  results  of  DIVPRK  are 
returned  to  the  subroutine  TEMP  as  a  temperature  time  history. 
TEMP  then  calculates  a  sum  of  the  squares  error  between  the 
experimental  data  read  into  the  program  PID  and  the  results  of 
the  equation  solver.  This  error  is  returned  to  the  optimizer 
program  where  the  optimizer  adjusts  the  unknown  parameters  and 
repeats  the  process  until  a  set  of  convergence  criteria  is 
met. 

Once  the  program  converges  the  results  are  vrritten  to  two 
files  Result, dat  and  Temp. mat.  Results.dat  contains  the 
calculated  values  of  the  convective  resistance  and  the 
convective  heat  transfer  coefficient.  The  temp. mat  file 
contains  the  experimental  and  model  temperature  time 
histories.  This  file  can  be  read  into  MATLAB  to  be  graphed, 
which  gives  a  visual  representation  of  how  close  the  match  is 
between  the  two  profiles. 

5.  Experimental  Data  Processing 

The  experimental  data  used  in  the  identification 
process  is  obtained  from  the  Naval  Weapons  Center  (NWC) ,  China 
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Lake,  California.  data  is  received  on  a  o  'ter  readout 
sheet  for  the  different  thermocouple  locations  on  the 
e*<perimental  vane.  The  temperature  is  measured  in  degrees 
Fahrenheit .  The  computer  program  needs  the  temperatures  in 
degrees  kelvin  and  referenced  to  the  ambient  or  reference 
temperature.  This  referencing  of  the  temperature  profiles 
enables  the  initial  conditions  of  the  equation  solver  to  be 
set  to  zero.  In  order  to  use  the  experimental  data  it  must 
first  be  preprocessed.  The  preprocessing  of  the  data  is 
illustrated  in  the  block  diagram  in  Figure  3.3. 


Figure  3.3  Experimental  Data  Preprocessing  Diagram 
The  Pre.dat  file  is  first  produced  by  manually  selecting  the 
series  of  tenperatures  of  a  specified  time  increment.  The 
separation  of  the  temperatures  is  up  to  the  individual.  The 
time  interval  selected  when  recording  the  experimental 
temperatures  must  also  be  placed  into  the  TEMP  subprogram  for 
the  differential  Equation  solver  DIVPRK.  This  time  interval 
allows  the  eciuation  solver  to  calculate  the  temperature  time 
histories  over  the  S6une  time  increments  as  the  experimental 
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data  being  used.  Failure  to  match  these  time  increments  may 
cause  erroneous  values  when  trying  to  match  experimental  and 
calculated  values.  The  values  are  entered  into  the  file 
Pre.dat  manually  as  a  nxl  matrix.  The  raw  data  from  the 
Pre.dat  file  is  read  into  a  preprocessing  program  called 
Trans. for.  This  smoothing  program  involves  adding  the  two 
ten^eratures  prior  to  the  one  being  averaged,  the  temperature 
being  averaged  and  the  two  temperatures  subsequent  to  that 
particular  tenperature.  Since  there  are  no  temperature  values 
prior  to  the  first  two  and  subsequent  to  the  last  two,  these 
values  are  taken  as  is.  For  61  temperature  values  over  a 
given  period  of  time  temperatures  one  and  two,  and  60  and  61 
are  taken  as  is.  The  temperatures  3  through  59  are  averaged 
values  once  the  Trans. for  program  has  been  run. 

The  number  of  files  produced  from  Trans. for  may  vary 
depending  on  the  number  of  experimental  temperature  profiles 
used.  The  basic  files  produced  are  Datam.dat  and  TV. mat.  The 
IM.mat  file  is  only  used  if  experimental  data  is  being 
utilized  from  the  mount.  The  files  TV. Mat  and  IM.mat  are  read 
into  a  MATLAB  program  called  Pro.m.  This  program  normalizes 
the  tenqperatures  to  the  ambient  and  converts  the  temperature 
from  degrees  fahrenheit  to  degress  Kelvin.  The  output  of  this 
program  is  a  file  called  Datam.dat  which  is  a  nxl  matrix. 
Datam.dat  is  the  smoothed  file  containing  the  time*  temperature 
profiles  of  the  experimental  data.  At  this  point  the 
Datam.dat  file  is  ready  to  be  used  by  the  main  program. 
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6.  Forward  Model  Matching 

a.  Coapsir±son  of  Steady  State  T&igperatiires 

A  forward  model  is  the  reverse  of  the  vane 
optimization  model.  It  takes  a  set  of  values  and  generates  a 
set  of  time-temperature  profiles  for  a  given  set  of  data. 
There  is  no  optimizing  of  the  variables.  All  the  rate 
coefficients  are  constants  and  the  profiles  generated  are 
based  on  the  solution  to  the  differential  equations  using 
these  constants  over  a  given  period  of  time.  A  forward  model 
was  constructed  from  the  program  PID  to  check  the  types  of 
time-tenperature  profiles  that  would  be  generated  by  these 
equations  given  a  set  of  values  for  the  rate  coefficients.  To 
construct  the  forward  model  from  the  main  program  PID  the 
optimizer  subroutine  was  eliminated  and  the  PID  program  called 
the  TEMP  subroutine  directly.  This  forward  model  allows  the 
user  to  solve  a  set  of  equations  over  a  given  period  with  a 
set  of  rate  coefficients,  and  b^^,  remaining  constant.  A 
set  of  time-tenperature  profiles  was  constructed  from  the  full 
scale  energy  balance  equations  using  the  final  optimized 
values  for  the  unknowns  that  Reno  found  from  the  quarter  scale 
modeling  process  [Ref.  l:p.  10].  The  values  were  scaled  up 
from  quarter  to  full  scale  and  inserted  into  the  PID  program 
for  the  unknown  values  at  the  corresponding  locations.  The 
values  of  b22  and  b24  did  not  correspond  to  anything  that  was 
done  by  Reno  (Ref.  1],  so  the  values  for  the  convective  rate 
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coefficients  were  all  set  equal  to  each  other.  The  forward 
model  was  run  on  the  six  node  vane  model  for  both  quarter  and 
full  scale  using  the  following  values:  Quarter  Scale  - 
a34=0.299,  a43=0.1322,  a4c=0.1018,  bu=30.163,  b22=0.271  Full 
Scale  -  a56=0. 01869,  a65=0. 008262,  aeG=0 .006362,  bii=1.885,  b22= 
b32=  b42  =0.017.  If  the  thrust  profile  is  allowed  to  continue 
indefinitely  the  steady  state  values  of  tenperature  should  be 
the  same  for  both  the  quarter  and  full  scale  values  on  the  six 
node  vane  model.  These  values  can  be  determined  analytically 
by  considering  a  voltage  divider  across  the  resistance  path  to 
ground.  The  convective  temperature  sources  act  as  the  supply 
voltages  in  such  a  case.  The  analytical  calculations  for  both 
sets  come  up  to  the  same  steady  state  temperature  values  for 
both  scaling  sizes  as  shown  in  Figure  3.4.  The  main 
difference  is  that  the  quarter  scale  reached  its  steady  state 
ten^erature  at  a  much  quicker  rate  than  the  full  scale.  The 
quarter  scale  reached  steady  state  in  approximately  20  seconds 
where  as  the  full  scale  had  not  reached  its  steady  state 
ten:«)eratures  even  after  160  seconds.  Figure  3.4  shows  the 
quarter  scale  response  of  the  node  temperatures  over  a  160 
seconds  firing  period. 
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Figure  3.4  Quarter  Scale  Forward  Model  Tai^erature  Response 
The  response  curves  in  Figure  3.4  are  for  a  six  node  vane 
model  with  the  nodes  at  the  shaft  and  mount  being  neglected. 
The  analytical  values  of  the  steady  state  temperatures  at 
nodes  one  through  four  are;  Tl=2340  K.  T2a2210  K,  T3al639  K 
and  T4s926  K,  obtained  from  the  figure  above.  The  values  of 
the  rate  coefficients  are  set  to  those  corresponding  to  the 
full  scale  vane  and  the  forward  model  program  is  executed 
again.  The  time  temperature  profiles  for  the  full  scale  are 
listed  in  Figure  3.5  over  the  same  160  second  firing  period. 
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Figure  3.5  Full  Scale  Forward  Model  Tesoperature  Response 
The  analytical  calculations  yield  the  same  steady  state  values 
for  nodes  one  through  four.  The  variation  is  in  the  amount  of 
time  needed  to  get  to  the  final  value.  The  full  scale 
requires  a  much  greater  period  of  time  to  reach  the  steady 
state  than  the  quarter  scale  does,  although  the  full  scale 
will  reach  the  same  steady  state  temperature.  This  goes  back 
to  the  time  constants  of  the  two  different  scale  models.  The 
rate  characteristics  decrease  for  the  full  scale  model  but 
these  are  the  reciprocals  of  the  time  constant  which  actually 
increase  for  the  full  scale  model. 
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b.  Rasponse  Time  Characteristics 

A  check  needs  to  be  made  to  see  if  the  response 
characteristics  of  the  energy  balance  equations  for  the  six 
node  model  are  similar  to  those  of  the  four  node  model 
obtained  from  Reno's  work  [Ref.  1] .  The  forward  model  will  be 
utilized  with  the  values  of  the  rate  coefficients  set  equal  to 
those  in  Reno's  thesis  [Ref.  l:p.  10].  The  thrust  profile 
will  be  the  same  as  that  used  by  Reno  and  so  will  the  recovery' 
temperatures,  with  Tri=2360  and  Tr2  =2260.  The  time- temperature 
profiles  obtained  by  Reno  for  a  four  node  quarter  scale  model 
are  shown  in  Figure  3.6. 
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Figure  3.6  Temperature  Histories,  Four  Mode  Model 
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The  results  in  Figure  3.6  are  based  on  an  optimization  program 
but  if  the  energy  balance  equations  are  correct  the  final 
values  used  by  Reno  should  be  able  to  be  inserted  into  the 
forward  model  generating  a  set  of  temperature  profiles  similar 
to  those  shown  above.  Figure  3.7  shows  the  results  of  this 
forward  model  using  the  six  node  energy  equations  with  Reno's 
rate  characteristic  coefficients. 


TUIE  (••€) 

Figttre  3.7  Six  Node  Quarter  Scale  Forward  Model  Response 
Comparison  of  the  temperature  profiles  of  Figures  3.6  and  3.7 
indicate  they  are  very  similar  with  the  temperature  profiles 
of  T2  from  the  four  node  and  T3  from  the  six  node 
corresponding  to  each  other.  T3  from  the  four  node 
corresponds  to  T5  from  the  six  node  and  T4  from  the  four  node 
corresponds  to  T6  of  the  six  node  model . 
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7.  Model  Simulation  ZZSSQ  Optimizer 

a.  UnknocmB  and  ABBvapticma 

The  first  series  of  runs  utilized  an  optimizer 
knovm  as  ZXSSQ  which  is  an  old  IMSL  routine  that  is  stored  in 
a  separate  library.  It  must  be  noted  at  this  point  that  the 
thrust  profile  for  the  full  scale  modeling  process  is  not  a 
ramp  function  as  those  conducted  in  the  quarter  scale  process 
but  is  being  simulated  as  a  step  response.  The  thrust  turns 
on  at  approximately  time  zero  and  shuts  off  approximately  6.7 
seconds  later.  Prior  to  starting  the  PSI  process  using  the 
program  formulation  provided,  the  unknowns  must  be  identified 
from  the  energy  balance  equations  and  a  set  of  assumptions 
must  be  established.  The  unknown  values  that  need  to  be 
determined  are  the  convective  resistances  from  the  free  stream 
to  the  vane,  Rp22«  Rp23'  %34,  conductive  resistance 
between  nodes  five  and  six,  Rgg,  the  conductive  resistance 
from  node  six  to  ground,  R^q,  and  the  capacitance  of  node  six, 
Cg.  Although  the  tip  convective  resistance  and  bj^j^  are 
unknown,  the  first  series  of  the  identification  process 
assumed  this  value  to  be  constant  as  it  was  in  the  quarter 
scale  siimilation  conducted  by  Driels  [Ref.  3}  and  Reno  (Ref. 
11. 

Two  simulations  were  conducted;  the  first  with  the 
assumption  that  the  convective  resistances  to  the  free  stream 
were  dependent  and  equal,  that  is  Rp22*%23"%24'  makes 
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b22=b32=b42*  The  second  simulation  made  the  assumption  that  the 
convective  resistances  were  independent  and  not  equal. 

Jb.  Utilization  of  Experimantal  Vane  Tea^  Profiles 

Prior  to  model  simulation  the  node  locations  for 
the  use  of  experimental  data  needs  to  be  selected.  The 
quarter  scale  modeling  done  in  previous  work  has  utilized 
experimental  ten^erature  profiles  from  the  shaft  and  mount 
locations.  For  the  full  scale  vane  experimental  data  is 
available  from  the  fin  area  of  the  vane.  A  set  of  simulations 
will  be  conducted  using  the  same  set  of  locations  as  that  done 
by  the  quarter  scale  modeling  process.  Then  the  experimental 
data  will  be  changed  and  the  temperature  profiles  from  nodes 
two,  three  and  four  will  be  utilized  as  the  experimental  time- 
ten^erature  profiles  to  be  matched. 

c.  Results 

A  series  of  simulations  were  conducted,  with  the 
initial  conditions  set  at  zero,  varying  the  assumption  of 
dependence  and  independence  of  the  convective  rate 
coefficients  and  varying  the  experimental  data  used  in  the 
identification  process.  All  the  simulations  had  the  same 
outcome,  there  was  no  convergence.  The  values  of  as«  and 
increase  substantially  in  value  causing  a  fatal  error  three  in 
the  differential  equation  solver  giving  the  error  message  that 
the  differential  equation  problems  became  too  stiff.  An 
attempt  was  made  to  alter  the  starting  initial  conditions  in 
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the  hope  that  the  values  could  be  placed  closer  to  their  final 
values.  This  was  a  manual  attenpt  to  correct  the  sudden 
change  in  the  values  for  the  unknowns.  Varying  the  initial 
conditions  had  no  effect  on  the  outcome.  The  problem  still 
failed  to  converge. 

8.  Model  Simulation  DBCLSF  Optimizer 
a.  UsB  of  Constrained  Optimizer 

Since  the  values  of  some  of  the  unknowns  would  jump 
substantially  causing  the  program  to  fail,  a  constrained 
optimizer  was  used  in  an  attempt  to  hold  these  values  within 
a  certain  set  of  limits.  The  constrained  optimizer  chosen  was 
DBCLSF  which  is  a  current  IMSL  subroutine.  The  reasoning 
behind  this  was  that  maybe  the  jump  in  the  unknowns  was  caused 
by  a  local  minima  and  when  the  optimizer,  ZXSSQ,  made  a  large 
adjustment  it  caused  the  problem  to  become  to  stiff.  The 
constrained  optimizer  would  keep  the  values  of  the  unknowns 
low  enough  to  prevent  the  differential  equation  solver  from 
becoming  to  stiff  and  allow  it  to  continue  to  solve  the  set  of 
equations  given. 

Jb.  Boundary  Limitations 

The  constrained  optimizer  allows  the  user  to  select 
a  variety  of  boundary  conditions.  This  optimizer,  DBCLSF, 
permits  the  selection  of  boundary  limits  that  are  all 
negative,  all  positive,  or  manually  set  either  to  each 
individual  unknown  or  as  a  blanket  setting  for  all  the 
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unknowns.  A  simulation  was  conducted  with  the  boundary  limits 
set  to  all  positive  and  then  to  a  set  with  the  lower  boundary 
being  set  at  0.0001  and  the  upper  boundary  set  at  100. 

c.  Rasulta 

For  the  simulation,  with  the  set  of  initial 
conditions  at  zero,  there  was  no  convergence  with  either  the 
boundary  limits  set  all  positive  or  for  the  upper  and  lower 
boundary  limits.  Again  the  initial  starting  values  of  the 
unknowns  were  changed  from  zero  to  various  values  in  an 
attempt  to  avoid  any  possible  local  minima.  The  results  were 
the  same  as  before  the  program  failed  to  converge.  The 
program  had  an  arithmetic  fault  or  just  continued  on 
indefinitely  with  the  values  of  the  unknowns  reaching  the 
boundary  limits  and  remaining  there.  The  program  is  trying  to 
match  a  series  of  time-temperature  profiles  from  the  fin  area 
of  the  vane.  The  program  does  not  converge  due  to  the 
substantial  jun^)  in  the  unknowns  at  the  mount  location.  From 
the  experimental  data  the  temperature  profiles  at  the  shaft 
and  mount  locations  vary  little  from  the  ambient.  Elimination 
of  the  mount  node,  node  six,  should  have  little  effect  on  the 
predictions  for  nodes  two,  three  and  four.  This  coupled  with 
the  inconsistencies  in  use  of  mount  configurations  leads  to 
the  modification  of  the  original  six  node  model  to  a  five  node 
model.  A  model  without  a  mount  node  and  with  the  shaft  node, 
node  five,  going  directly  to  ground. 
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B.  FIVE  MODE  MODEL  APPROACH 

1.  Concept  Behind  Mbdlfication  From  Six  Node  Model 

The  computer  simulations  conducted  on  the  six  node 
model  encountered  the  same  problem  in  each  simulation.  The 
values  of  the  unknowns  ass  and  a^s,  corresponding  to  the  mount 
location,  would  juir^)  a  large  amount  in  value  causing  the 
program  differential  equation  solver  to  experience  a  fatal 
error,  indicating  that  the  equations  were  too  stiff.  Changing 
the  relationship  between  the  convective  rate  coefficients,  bj^, 
had  no  effect  on  the  resultant  outcome  and  the  jump  in  these 
unknowns  prevented  the  convergence  of  the  program.  These 
unknowns  contribute  to  the  tenperature  profiles  primarily  at 
the  mount  location.  The  convective  heat  transfer  coefficients 
that  were  to  be  identified  were  between  the  free  stream  and 
nodes  two,  three  and  four  of  the  vane  and  at  this  point  we 
were  not  concerned  with  the  mount.  Examination  of  the 
experimental  time  tenperature  profiles  at  the  mount  location, 
for  the  full  scale  vane,  revealed  that  the  temperature  at 
these  locations  had  little  variation.  The  temperature  varied 
less  than  25°F  over  the  entire  >.7  second  rocket  firing  period 
as  conpared  to  the  temperature  variation  at  nodes  two,  three, 
and  four  which  varied  from  900  to  3000®  F  during  the  same 
rocket  firing  time.  In  contrast  to  the  quarter  scale  modeling 
the  temperature  variation  in  the  mount,  node  four,  was 
significant.  The  size  of  the  quarter  scale  vane  limited  the 
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locations  for  attachment  of  thermocouples  to  the  shaft  and  the 
mount.  This  limited  the  experimental  temperature  data  to  only- 
two  locations,  the  mount  and  shaft,  making  it  necessary  to  use 
the  mount  temperature  profiles  in  the  matching  process.  The 
full  scale  vane  is  larger  in  size  and  provides  greater  space 
to  place  thermocouples  and  collect  the  experimental  data. 

The  model  was  constructed  to  determine  the  convective  rate 
coefficients  from  the  free  stream  to  the  vane.  Due  to  the 
small  variation  of  temperature  in  the  mount  and  the  enormous 
change  of  the  unknowns  using  the  mount  node  it  was  decided 
that  a  modification  needed  to  be  made  to  the  six  node  model. 
This  modification  was  the  elimination  of  the  mount  node,  which 
eliminated  the  two  unknowns,  age  and  a^s  that  were  causing  the 
program  to  fail.  There  was  also  another  motivation  for 
eliminating  the  mount  node,  which  was  the  inconsistency  in  the 
use  of  the  mount  apparatus  for  all  the  experimental  rocket 
firings  (Ref.  l:p.  13] .  The  modification  to  a  five  node  model 
is  illustrated  in  Figure  3.8. 
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Figure  3.8  Five  Node  Vane  Model 
2.  Energy  Balance  Equation  Formulation 

The  formulation  of  the  energy  balance  equations  for 
the  five  node  model  was  done  using  the  same  process  as  that 
for  the  six  node  with  only  minor  changes.  The  modification 
involves  eliminating  node  six  which  in  turn  eliminates  the 
last  equation  from  the  six  node  model  energy  balance.  There 
was  also  only  a  minor  change  to  the  equation  for  node  five  of 
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the  original  six  node  equations.  The  following  equations  are 
applicable  to  the  five  node  vane  model; 


i*! = Ty +ai2  ^2  Tgj 

(45) 

^2  ~^2l  ~^22  ^2  ■‘■^23  ^3  '*■■^2 

(46) 

^3  "^32^2  “^33  ^3  '*■^34  ^4  ■*‘■^2 

(47) 

i’g  =843  Tj  -844  +845  Tg  +b42  Tg2 

(48) 

^5“®54^4~®55^5 

(49) 

C^R^g  ^55“^54'‘‘350 

(50) 

Reduction  of  the  model  to  five  nodes  eliminates  the  need  to 
identify  the  unknown  conductive  resistances  from  the  shaft  to 
the  mount,  the  mount  to  ground,  Rgo,  and  the  capacitance 
of  the  mount  node,  Cg.  The  experimental  data  for  the  five 
node  full  scale  modeling  process  will  be  from  the  thermocouple 
attachment  points  that  correspond  to  nodes  two,  three  and  four 
of  the  vane  model. 

3.  Program  Modification 

The  main  program  PID  needs  to  be  modified  from  the  six 
node  model  to  that  of  the  five  node.  This  is  done  primarily 
by  making  an  alteration  in  the  FCN  subroutine  which  generates 
the  differential  equations  to  be  used.  The  only  other 
modification  that  needs  to  be  made  to  the  PID  program  is  the 
set  up  of  the  size  of  the  dimensions  for  the  rate  coefficients 
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and  the  setting  of  the  initial  conditions.  Identification  of 
the  unknowns  will  be  discussed  in  the  following  section.  The 
FORTRAN  coding  used  to  generate  the  simulation  program  was 
constructed  to  be  flexible  in  that  it  allows  for  changes  in 
the  vane  nodes  without  major  program  modification. 

4.  Model  Simulation  ZXSSQ  Optimizer 
a.  Uziknovns  and  Asstaaptxons 

The  modeling  process  for  the  five  node  model  will 
be  set  up  in  the  same  way  the  six  node  model  utilized  the 
ZXSSQ  optimizer.  Prior  to  starting  the  simulations  the 
unknowns  for  this  model  need  to  be  identified.  As  before  the 
unkiiowns  are  the  convective  rate  coefficients,  b22,  b32  and  b42, 
along  with  the  conductive  rate  coefficient,  In  the 
quarcer  scale  modeling  it  was  demonstrated  that  the  convective 
rate  coefficient  from  the  tip,  b^,  had  little  effect  on  the 
determination  of  the  unknown  parameters.  Without  knowing  the 
effect  on  the  full  scale  model  the  value  of  b^i  will  be 
considered  a  constant  for  the  first  series  of  simulations. 
The  assumption  of  dependence  between  the  convective  rate 
coefficients  was  applied  to  the  first  series  of  simulations 
making  b33=b}2=b42.  A  second  set  of  simulations  will  be 
conducted  making  the  values  of  the  convective  rate 
coefficients  independent  of  each  other  and  making  the  tip 
convective  rate  coefficient  a  variable. 
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b.  Ea^riaental  Data  Used 

The  experimental  data  used  for  the  five  node  vane 
model  will  be  from  the  thermocouple  locations  that  correspond 
to  nodes  two,  three,  and  four  of  the  model.  For  the 

simulations  involving  the  use  of  the  ZXSSQ  optimizer,  an 
attempt  will  be  made  to  match  all  three  time-tenperature 
profiles.  A  point  to  be  noted  is  that  the  values  of  the 
recovery  temperatures  that  are  to  be  utilized  are  different 
from  those  used  in  the  quarter  scale  model  due  to  the 
different  rocket  motor  used  in  the  test  firings.  The 
tenperatures  used  to  drive  the  heat  transfer  process  were 
Tri=2970  and  Tr2=2870  Kelvin,  respectively  [Ref.  2:p.  44]. 
Referencing  these  tenperatures  to  the  ambient,  as  was  done 
with  the  experimental  node  tenperatures,  would  lead  to 
recovery  temperatures  being  normalized  to  values  of  Tri=2670 
and  Trj»2570  Kelvin, 
o.  Results 

As  seen  with  the  six  node  modelling  process,  using 
the  ZXSSQ  optimizer,  the  model  failed  to  converge.  Changing 
the  assumption  for  the  dependence  of  the  convective  rate 
coefficients,  which  made  the  values  of  the  convective  rate 
coefficients  independent  variables  had  no  effect,  and  the 
program  still  failed  to  converge.  An  attenpt  was  made  to 
match  a  single  node  time  temperature  profile  but  this  also 
proved  fruitless  and  again  the  program  failed  to  converge. 
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From  the  foregoing  results  it  is  evident  that  the  use  of  the 
ZXSSQ  optimzer  for  the  full  scale  modeling  process  is 
inadequate.  There  is  no  provision  in  the  ZXSSQ  optimizer  to 
control  the  values  or  range  of  the  unknown  variables. 

5.  Mbdel  Simulation  DBCLSF  Optimizer 
a.  Conatra.in3d  Modal  Simulation 

As  with  the  six  node  model  simulations,  a 
constrained  optimizer  replaced  the  present  optimizer,  ZXSSQ. 
This  was  an  attempt  to  hold  the  unknowns  within  certain  limits 
permitting  the  optmizer  a  chance  to  avoid  or  recover  from  any 
possible  local  minimum  and  prevent  the  unknowns  from  becoming 
so  large  that  the  program  fails  to  converge.  The  DBCLSF 
constrained  optmizer  utilized  will  be  from  the  IMSL  library. 
The  first  simulation  assumed  that  the  convective  rate 
coefficients  bjj,  Dj,  and  b„  are  equal  and  that  b^  was  a 
constant  and  attempted  to  match  the  ten^erature  profiles  of 
nodes  two,  three,  and  four.  This  simulation  also  failed  to 
converge.  The  assumption  that  bn  was  a  constant  was  changed 
and  bxi  was  allowed  to  vary.  This  simulation  also  failed  to 
converge . 

b*  Single  Soda  Temperature  Pxofila  Modaliag 

A  different  approach  was  taken  to  try  and  establish 
a  five  node  working  model  of  the  vane.  Instead  of  trying  to 
match  all  three  temperature  profiles  the  model  was  altered  to 
try  and  match  only  a  single  node  time-teitq;>erature  history  then 
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work  up  to  two  nodes  and  then  a  full  three  node  matching 
process.  The  assunptions,  for  the  single  node  tenperature 
matching  process,  that  any  of  the  unknown  variables  were 
dependent  were  abandoned  and  each  unknown  variable  was 
considered  independent  of  the  others. 

(1)  Node  Two  Tenperature  Modeling.  The  modeling 
simulations  began  by  trying  to  match  the  time-tenperature 
history  of  only  a  single  node,  node  two.  This  was  attenpted 
in  the  hope  that  if  only  one  profile  could  be  matched  that  by 
examination  of  the  profiles  generated  for  the  other  nodes  a 
slight  adjustment  could  be  made  in  order  to  obtain  a  match 
between  all  three  node  time  tenperature  profiles.  The  program 
was  modified  to  match  the  profile  of  node  two  and  the  boundary 
limits  were  set  to  0.0001  for  the  lower  boundary  and  100  for 
the  upper  boundary,  with  recovery  temperatures  of  T*i»2670  and 
Tm»2570,  and  the  initial  conditions  set  to  zero.  The 
simulation  converged  with  a  sum  of  the  squares  error  of  12.01, 
with  values  for  the  unknowns  b„=3.0833,  b,, *0.0255,  b,a*0.0492, 
and  b^, *0.2496,  The  value  for  a«  reached  the  lower  limit  of 
0.0001  set  in  the  constrained  optimizer  and  did  not  vary  from 
that  point.  Figure  3.9  illustrates  the  time-temperature 
profiles  for  this  sisoilation. 
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Figure  3.9  Time-Tenperature  Profiles  Node  Two  Match 

(2)  Node  Three  Temperature  Modeling.  The  main 
program  PID  was  modified  to  match  experimental  data  for  node 
three  only.  The  program  simulation  converged  with  a  sum  of 
the  squares  error  of  15.26.  The  values  obtained  for  the 
variables  were  bu=2.7883,  b225=0.8934,  bjjsO.OOOl  and  b42=0 . 1778 . 
Again  the  value  for  as<,  dropped  to  the  lower  limit  allowed  by 
the  constrained  optimizer,  0.0001,  and  remained  there.  Figure 
3.10  illustrates  the  time- temperature  profiler  of  the  node 
three  match. 
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T-gure  3.10  Time-Tesnperature  Profiles  Node  Three  Match 

(3)  Node  Four  Temperature  Modeling,  Again  the 
program  was  modified  to  match  the  experimental  data  for  node 
four  only  and  the  simulation  converged  for  this  set  up  with  a 
sum  of  the  squares  error  of  39.539.  The  values  obtained  for 
the  unknowns  were  bus2.92,  bj2s»1.38,  bj, *0.212  and  b4i=0.02  and 
again  the  value  of  a^c  went  to  the  lower  boundary  limit  sec  in 
the  constrained  optimizer  and  remained  there.  The  time- 
temperature  profiles  are  illustrated  in  Figure  3.11. 
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Figure  3.11  Time-Tenperature  Profiles  Node  Four  Match 

(4)  Results.  In  all  three  simulations,  using  only 
a  single  set  of  experimental  data,  the  program  converged  and 
the  calculated  and  experimental  profiles  of  the  node  that  was 
being  matched  were  as  accurate  as  the  optmizer  could  get.  The 
error  using  node  two  data  was  less  than  that  of  node  three  and 
much  less  than  that  of  node  four.  The  problem  arises  in  the 
calculated  temperature  histories  of  the  other  two  nodes  not 
being  matched.  It  can  be  seen  in  Figures  3.9,  3.10  and  3.11 
that  the  two  nodes  not  being  matched  are  significantly 
different  from  the  experimental  data.  In  the  simulation  using 
the  experimental  data  from  node  two  the  profiles  of  node  two 
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and  three  match  fairly  closely  while  the  profile  for  node  four 
is  nowhere  near  the  experimental  data.  In  the  simulation 
using  the  node  three  experimental  data,  node  three  profiles 
have  a  fair  match  but  the  profiles  for  nodes  two  and  four  are 
unsatisfactory.  The  same  condition  holds  when  using  the 
experimental  data  for  node  four  only.  It  must  be  noted  that 
the  magnitude  of  the  error  between  the  calculated  and 
experimental  data  is  dependent  on  the  quality  of  the 
experimental  data.  If  the  data  follows  a  first  order  linear 
response  closely  then  the  match  of  the  calculated  and 
experimental  profiles  is  adequate  and  the  sum  of  the  squares 
error  is  low.  If  the  experimental  profiles  vary  from  an 
exponential  curve  then  the  first  order  linear  model  matches 
the  experimental  profile  the  best  it  can  and  the  error  will  be 
large . 

o.  Muiual  Adjustmmnt  of  Vmriablma-Formrd  Hodal 
The  results  of  the  single  node  tenperature  profile 
simulation  permit  the  use  of  the  forward  model  to  generate  a 
series  of  time- temperature  profiles  with  the  final  values 
obtained  for  the  unknowns.  By  adjusting  the  various  unknowns 
we  can  see  the  effect  one  variable  has  on  another.  A  set  of 
values  needs  to  be  selected  for  the  unknowns.  These  values 
can  be  arbitrary  since  we  are  looking  at  the  effects  of  one 
value  on  another.  The  initial  starting  values  for  the 
unknowns  were  chosen  as,  busi.ve,  b2as0.0455,  bj2=>0.06,  b4i=0.02 
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and  aso-O.OOOl.  The  value  of  aso  was  chosen  as  0.0001  because 
in  all  the  previous  simulations  this  unknown  juitped  to  this 
lower  value  limit  and  remained  there.  The  values  for  the 
other  unknowns  were  chosen  by  scaling  up  the  final  values  of 
the  quarter  scale  model  results  obtained  by  Reno  [Ref.  1]  and 
considering  the  values  obtained  from  the  single  node 
teirperature  profile  results.  The  effects  of  the  variables 
were  that,  the  value  of  bu  influenced  the  teirperature  profile 
of  node  two  and  had  little  or  no  effect  on  nodes  three  and 
four,  while  the  same  holds  true  for  the  other  convective  rate 
coefficients.  Altering  each  one  consecutively  varied  the 
result  of  that  particular  node  profile  but  had  little  effect 
on  the  outcome  of  the  other  node  calculated  temperature 
profiles.  As  an  example  of  this  effect  the  values  of  the 
unknowns  listed  above  were  used  in  the  forward  model.  The 
model  used  a  set  of  recovery  tenperatures  with  values  of 
Tri=2970  and  Trj=2870.  There  was  no  need  in  normalizing  these 
recovery  tenperatures  in  the  forward  model.  The  effects  of 
one  variable  on  the  other  was  being  examined  and  not  actually 
trying  to  determine  the  proper  values  for  the  unknowns  as  done 
in  the  PID  modeling  program.  The  profiles  generated  by  this 
forward  model  simulation  are  illustrated  in  Figure  3.12  where 
the  three  nodes  of  concern,  nodes  two,  three,  and  four  are 
displayed . 


52 


“0  05  1  U  2  15  .  3  4.  5 

TIME  (sec) 

Figure  3.12  Forward  Model  Time-Tenperature  Profiles 
These  profiles  show  a  good  match  between  the  calculated  and 
experimental  ten^erature  profiles  of  node  two.  The  profiles 
of  nodes  three  and  four  do  not  follow  the  experimental  data 
but  again  the  experimental  data  for  these  two  nodes  does  not 
follow  an  exponential  curve.  Another  thing  that  needs  to  be 
pointed  out  is  the  calculated  curves  for  the  three  nodes 
appear  to  be  linear  straight  lines.  This  is  not  an 
exponential  curve  as  esqpected  from  a  first  order  linear  model. 
To  show  the  effects  of  changing  one  variable  on  the  other 
temperature  profiles,  the  value  of  b)}  was  changed  to  0.03  and 
the  forward  model  simulation  was  executed.  The  resulting 
time^temperature  profiles  are  displayed  in  Figure  3.13. 
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Figure  3.13  Forward  Xodel  Simula tlon  Node  3  Variant 
Comparison  o£  Figures  3.12  and  3.13  show  that  there  was  a 
significant  change  in  the  calculated  tenperature  profile  of 
node  three,  while  the  profiles  of  nodes  two  and  four  changed 
very  little  in  magnitude  cou^red  to  node  three.  The  program 
calculated  the  sum  of  the  squares  error  for  both  simulations. 
The  error  for  the  first  simulation  was  74.4  while  the  error 
for  the  second  simulation  was  213.577. 

Adjustments  can  be  made  to  each  of  the  unknown 
variables  in  the  forward  model  to  create  a  better  match  for 
any  of  the  Individual  curves.  This  process  of  adjusting  the 
unknowns  is  a  manual  atten^t  to  do  what  the  optimizer  was 
intended  to  do  in  the  PIO  program.  The  manual  adjustment 
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allows  the  user  to  view  the  effects  that  individual  rate 
coefficients  have  on  the  temperature  profiles  of  the  other 
nodes . 

d.  Ttio  Node  T&uperature  Profile  Modeling 

The  previous  section  established  the  convergence  of 
the  modeling  program  and  the  matching  of  a  single  time- 
temperature  profile.  The  model  will  be  upgraded  to  matching 
two  temperature  profiles  in  an  atten^t  to  obtain  a  three  node 
full  scale  matching  model. 

(1)  Nodes  Two  and  Three  Temperature  Modeling.  The 
modeling  program  was  modified  to  match  the  temperature 
profiles  of  two  nodes,  vice  one,  with  the  sum  of  the  squares 
error  being  calculated  based  on  both  temperature  profiles. 
The  first  simulation  will  attenpt  to  match  the  temperature 
profiles  of  nodes  two  and  three.  The  assumptions  for  the  two 
node  matching  simulations  are  the  same  as  those  for  the  single 
node  in  that  all  the  unknowns  are  considered  independent  of 
each  other  and  the  constrained  optimizer  will  be  utilized  with 
the  upper  boundary  of  100  and  a  lower  boundary  of  0.0001.  The 
simulation  converged  with  the  resulting  sum  of  the  squares 
error  of  15.217  and  values  for  the  unknowns  bii-2.83, 
b22**0.0297,  b32'"0.0056,  b42'<‘0.9408  and  asg-O.OOOl.  Again  the 
value  for  a^Q  jumped  to  the  lower  limit  and  remained  there 
while  the  optimizer  was  changing  the  values  for  the  other 
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unknowns.  The  time -temperature  profiles  for  the  two  node 
matching  process  are  illustrated  in  Figure  3.14. 


Figure  3,14  Two  Nod«  Tcoporaturo  Profllos 
The  simulation  program  matched  the  time-ten^erature  profiles 
for  nodes  two  and  three  while  the  calculated  profile  for  node 
four  was  no  where  near  the  experimental  profile.  Figure  3.14 
illustrates  that  the  calculated  profile  for  node  four,  T4 
calculated  profile,  was  unsatisfactory. 

(2)  Nodes  Two  and  Four  Temperature  Modeling.  The 
program  was  again  modified  to  match  the  data  for  the 
temperature  profiles  from  nodes  two  and  four.  The  simulation 
using  these  experimental  profiles  yielded  a  sum  of  the  squares 
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error  of  216.02  and  values  of  the  unknowns  1.2010, 
b22-0.1151,  b32-0.0001,  b42-0.0001  and  agQ-O.OOOl.  The  time 
temperature  profiles  for  this  simulation  are  illustrated  in 
Figure  3.15. 


Flgiure  3.15  Tiro  Node  Twaporaturo  Prof  lias 
The  experimental  and  calculated  profiles  for  node  two  matched 
adequately,  but  the  calculated  profile  for  node  four  is 
substantially  off  from  the  experimental  data.  The  value  of 
®5G  reached  the  lower  limit  set  by  the  optimizer  and 

remained  there. 

(3)  Nodes  Three  and  Four  Temperature  Modeling.  The 
final  two  node  matching  simulation  was  conducted  by  modifying 
the  main  program  and  calculating  the  sum  of  the  squares  error 
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based  on  the  temperature  profiles  from  nodes  three  and  four. 
The  simulation  yielded  values  of  bii=2.19,  b22«0.77,  b32»0.06, 
b42»0.0001  and  agg-O.OOOl  with  a  sum  of  the  squares  error  of 
185.95.  As  can  be  seen  in  Figure  3.16  the  program  adequately 
matched  the  profiles  for  nodes  three  and  four  but  there  was  a 
significant  difference  in  the  calculated  and  experimental  time 
temperature  profile  of  node  two. 


(4)  j?esults.  The  program  was  able  to  match  the 
ten^erature  profiles  from  node  two  and  to  a  lesser  extent  node 
three,  but  had  difficulty  in  matching  the  profiles  for  node 
four.  This  is  due  to  the  fact  that  the  experimental  profile 
for  node  four  is  not  an  exponential  curve.  The  first  order 
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linear  modeling  program  provides  an  descent  attempt  to  model 
the  experimental  data  and  provides  as  close  an  estimation 
possible.  The  use  of  the  node  four  experimental  data  causes 
a  significant  increase  in  the  sum  of  the  squares  error. 
Although  it  may  be  possible  to  get  an  adequate  match  for  two 
temperature  profiles  the  third  is  by  far  not  satisfactory, 
e.  Convergent  Simulation  Model 

(1)  Three  Node  Temperature  Profile  Matching.  It 
has  been  demonstrated  that  the  modeling  program  can  match  one 
and  two  temperature  profiles.  The  process  is  now  upgraded  to 
model  all  three  ten^erature  profiles  of  the  vane,  that  is  the 
modeling  of  nodes  two,  three,  and  four.  This  is  the  primary 
reason  behind  the  modeling  of  the  full  scale  vane. 

(2)  Unknowns  and  ABBunptions .  The  model  simulation 

for  the  three  node  ten^erature  matching  will  assume  that  the 
unknown  variables,  b^^^,  b22«  ^321  b42  and  a^Q,  are  independent 
of  each  other.  Previous  five  node  model  simulations 
considered  a  constant  in  the  same  manner  as  the  quarter 
scale  modeling  did.  Changing  the  value  of  in  the  quarter 

scale  model  had  little  effect  on  the  values  of  the  other 
imknowns  .  This  fact  has  not  been  established  in  the  full 
scale  modeling  process  which  was  the  reasoning  behind  nvaking 
bj^]^  a  variable.  The  initial  conditions  will  be  set  to  zero 
for  all  five  unknowns  and  the  values  for  the  recovery 
temperatures  are  T}y^<*i2670  and  Tr2*2570  Kelvin.  For  this 
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simulation  the  sxam  of  the  squares  error  were  be  calculated 
based  on  all  three  temperature  profiles  and  the  boundary 
limits  for  the  constrained  optimizer  were  set  to  0.0001  for 
the  lower  boundary  and  100  for  the  upper  boundary. 

(3)  Results.  The  simulation  was  executed  for  a 
time  period  of  zero  to  five  seconds.  The  simulation  converged 
with  a  sum  of  the  squares  error  of  58.5  and  the  values 
obtained  for  the  unknown  variables  were  bj^j^»1.3787,  b22»0.0938, 
b32"0.0862,  b42»0.0484  and  asQ-O.OOOl.  As  with  the  single  and 
two  node  temperature  matching  the  value  for  a^Q  dropped  to  the 
lower  boundary  set  in  the  optimizer.  The  low  value  obtained 
for  a5g  relates  to  the  negative  value  obtained  for  the  similar 
unknown  a4Q  of  the  quarter  scale  model  from  Oriels  [Ref.  3] 
simulation.  The  constrained  optimizer  prevented  the  value  for 
a^  frcxn  going  negative  in  the  five  node  model  simulation. 
The  time'ten^erature  profiles  generated  by  this  simulation  are 
illustrated  in  Figure  3.17. 
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Ftgur«  3.17  fiv*  Hod*  Van*  ICodsl  Taap^rature  Blstorlta 

The  match  between  the  profile  labeled  T2,  for  node  two,  was  an 
adequate  match  as  it  was  in  some  of  the  previous  simulations. 
The  match  between  the  experimental  and  calculated  profiles  for 
nodes  three  and  four  were  not  as  close  as  node  two.  Again 
this  is  due  to  the  quality  of  the  esqterimental  data  taken,  if 
the  data  does  not  follow  am  exponential  rise  such  as  a  first 
order  linear  system,  then  the  optimizer  is  less  accurate  in 
matching  the  experimental  temperature  profiles. 

The  simulation  above  was  conducted  using  the  recovery 
temperatures  of  Tri«2670  and  Tg2-2570  Kelvin.  These  values 
were  obtained  from  past  work  conducted  by  Nunn  (Ref.  l]  and 
Reno  (Ref.  2J .  Since  there  was  no  direct  way  to  determine 
these  values  for  the  full  scale  oxperiiuental  rocket  motor 
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firings,  due  to  the  unavailability  of  firing  data,  the  values 
of  the  recovery  temperatures  were  altered  to  illustrate  the 
effects  on  the  model  simulation.  The  first  modification  was 
to  change  the  values  of  the  recovery  temperatures  to  the  value 
of  the  chamber  temperature.  This  would  make  Tj^j^=2970  and 
Ti^2“2870  degrees  Kelvin.  This  simulation  converged  with  a  sum 
of  the  squares  error  of  56.5,  The  values  of  the  unknowns  were 
bii=1.4747,  b22=0.0716,  b32=0.0761,  b42=0.0430  and  a5G=0.0001. 
The  time- temperature  profiles  were  very  similar  to  those  in 
Figure  3.17.  Another  simulation  was  conducted  with  the  values 
of  Tjjj^=3550  and  Tjj2“3450  degrees  Kelvin.  This  simulation 
yielded  an  error  of  54.5  and  values  of  bj^j^=*2.3045,  b22"0*033, 
b32“0.0621,  b42"0.0354  and  aggaO.OOOl.  The  time -temperature 
profiles  for  this  simulation  are  again  very  similar  to  those 
obtained  in  Figure  3.17.  This  would  be  expected  due  to  such 
a  small  variation  in  the  sum  of  the  squares  error. 

f.  Simulation  With  various  Initial  Conditions 

To  test  the  robustness  of  the  simulation  program 
the  initial  conditions  were  varied.  Changing  the  initial 
conditions  changes  the  starting  point  for  the  model 
simulation.  Variation  of  the  starting  values  for  the  unknowns 
yielded  the  same  final  values.  The  program  converged  to  the 
same  point  and  the  same  profiles  when  the  initial  conditions 
were  changed  to  different  values. 
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g.  Calculations  of  Heat  Transfer  Coefficients 

The  values  obtained  for  the  unknovm  variables,  for 
the  three  node  temperature  matching  simulation  with  the 
recovery  tenperatures  of  Tj^3^=2670  and  Tj^=2570  Kelvin,  were 
b2^^®1.3787,  b22’*0*052®i  b22“0 •  ,  b^2*® •  and  agQ=0.0001« 

The  convecuive  heat  transfer  coefficients  can  now  be 
determined  but  the  values  for  the  resistance  must  first  be 
determined.  These  calculations  yield  Rjf3^=>0.1044  K/W, 
Rp22-0.3101  K/W,  Rp23”0.2530  K/W  and  Rp24-0.3605  K/W  with  the 
values  for  the  corresponding  convective  heat  transfer 
coefficients;  ht-22027.5  W/ro^*K,  W/m^-K,  h^3=l057 

W/m^*K,  and  hv4"594  W/m^«K.  Where  hj.  is  the  convective  heat 
transfer  coefficient  for  the  tip  and  1x^2 »  ^^3  Nr4 
convective  heat  transfer  coefficients  from  the  free  stream  to 
nodes  two,  three  and  four  of  the  vane,  respectively. 

The  simulation  involving  the  use  of  Tr3^«2970  and  Tr2“28?0 
Kelvin  ^/ielded  the  values  of  bj^;^-1.4747,  b22“0,07l6, 

b32-d.076l,  b42«0.0430  and  a^Q-O.OOOl.  These  values  yield 
resistances  of  Rpi-0.0976  K/W,  Rp22-0>4062  K/W,  Rp23-0.287  K/W, 
and  Rp24-0.4058  K/W  with  the  resulting  convective  heat 
transfer  coefficients  of  h^«2356l.3  W/m**K,  w/m^*K, 

1x^3-932.86  W/m^*K  and  hv4-527.74  W/m^*K.  The  last  simulation 
involved  changing  the  recovery  temperatures  to  "^Rl  -3550  and 
Tr2  -34*30  Kelvin.  The  results  of  this  simulation  yielded 
b33-2,  1045,  b22"0*033,  b32*0.C621,  b42»0.03S4  and  asg-0.0001. 
The  rc««i8tance  values  obtained  frcoa  these  results  are 
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Rfi=0.0624  K/W,  Rf22=0.8815  K/W,  Rf23=0.3512  K/W  and  Rf24=0.4928 
K/W.  Using  these  resistances  the  values  for  the  convective 
heat  transfer  coefficients  are  ht«36819  W/m^^K,  hy2=404-6 
W/m2*K,  h^3=761.5  W/m^^K  and  h^4=434.5  W/m^^K. 


64 


DISCUSSION  OF  RESULTS 


The  full  scale  six  noae  model  atteir^Jted  to  directly  scale 
up  the  four  node  quarter  scale  vane  model.  The  decision  to 
add  the  two  extra  nodes  was  necessary  in  order  to  model  the 
larger  vane  and  was  believed  to  more  accurately  describe  the 
heat  transfer  process  for  the  larger  amount  of  material 
contained  in  the  full  scale  vane.  The  first  series  of 
simulations  involved  the  use  of  the  2XSSQ  optmizer  to  adjust 
the  unknowns  in  order  to  obtain  a  match  between  the 
experimental  and  calculated  time-temperature  profiles.  The 
unknown  rate  coefficients  for  the  six  node  model  were 
identified  as  b,,,  b)3«  b^,  as(,  a^B  and  aso.  The  assumptions 
were  that  b^  was  a  constant,  as  it  was  in  the  quarter  scale 
simulations,  and  that  the  convective  rate  coefficients,  bjj, 
bj]  and  bt^,  were  equal.  These  series  of  simulations  failed  to 
converge  due  to  the  substantial  variation  in  the  unknown 
values,  Sbc  a^s  and  a«,  which  contain  the  resistances  between 
the  shaft  and  mount  and  the  mount  and  ground.  The  values  of 
a^s  and  ato  also  contain  the  unknown  value  of  the  mount  node 
capacitance.  During  the  simulation  these  variables  jumped 
substantially  in  value  causing  the  differential  equations  to 
become  stiff.  The  juiqp  in  the  variables  were  both  in  the 
negative  and  positive  direction.  The  experimental  data  for 
these  simulations  was  taken  from  the  shaft  and  mount  locations 
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in  the  same  manner  that  the  experimental  profiles  were  Ui;;ed  in 
the  quarter  scale  modeling  process.  When  the  simulation 
failed  to  converge  the  experimental  data  was  changed  to  that 
relating  to  nodes  two,  three  and  four  of  the  six  node  model. 
In  the  physical  test  firings,  for  the  full  scale  vane,  the 
tenperature  profiles  corresponding  to  the  shaft  and  mount 
locations  had  minimal  variation  from  the  ambient.  It  was 
believed  that  the  simulation  program  had  difficulty  adjusting 
to  the  small  changes  in  the  mount  and  shaft  temperature 
profiles  while  trying  to  make  large  adjustments  to  the 
unknowns  from  nodes  two,  three  and  four.  That  is,  over  a  time 
increment  the  tenperature  change  was  small  in  the  shaft  and 
mount  but  was  large  in  the  tip  and  fin  area.  This  idea  led  to 
abandoning  the  uf'  of  the  experimental  data  from  the  shaft  and 
mount  and  using  the  experimental  data  from  the  fin  area  nodes, 
nodes  two,  tnree  and  four.  The  simulation  using  this  new  set 
of  experimental  profiles  also  failed  to  converge  with  the  same 
results  for  the  unknowns,  a^^,  a^^  and  a(o> 

^wther  idea  considered  controlling  the  simulation  by 
preventing  the  unknown  values  from  varying  substantially. 
This  was  achieved  by  the  use  of  the  constrained  optimizer, 
which  would  permit  the  setting  of  boundary  limits  to  hold  the 
unknowns  within  a  certain  range.  This  range  was  selected  as 
positive  values  with  a  lower  boundary  limit  of  0.0001  and  an 
upper  boundary  limit  of  100.  The  assumptions  for  the  unknowns 
remained  the  same  as  the  previous  simulations.  This 
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simulation  using  the  constrained  optimizer  also  failed  to 
converge,  where  the  values  of  the  un^novms  asg,  a^s  and  agg 
reached  the  set  limits  and  remained  there.  Attenpts  were  made 
at  raising  the  upper  boundary  in  the  hope  that  the  higher 
boundary  limit  would  allow  the  optimizer  a  chance  to  avoid  or 
recover  from  any  local  minimum  and  continue  on  with  the 
simulation.  Raising  the  upper  boundary  limit  had  the  same 
result  as  before,  the  unknown  values  reached  the  boundary 
limits  and  remained  there  and  the  program  failed  to  converge. 
Any  attenpt  at  setting  the  lower  bou.ndary  limit  to  a  negative 
value  or  raising  the  upper  limit  to  too  high  a  value  caused  an 
arithmetic  fault  to  occur.  A  set  of  simulations  was  conducted 
using  each  optimizer  but  changing  the  assumptions,  making  the 
convective  rate  coefficients  b^],  bjj,  and  independent  of 
each  other.  In  both  cases  these  simulations  failed  to 
converge. 

Failure  to  resolve  the  problem  of  the  unknowns  associated 
with  the  shaft  and  mount  led  to  the  modification  of  the  model 
from  six  nodes  to  five.  This  was  done  by  eliminating  the 
mount  node  which  eliminated  the  unknowns  as«,  a«s  and  a^o  and 
added  the  unknown  ajo.  This  was  justified  by  the  fact  that 
there  was  no  consistent  mount  arrangement  used  from  one  test 
firing  to  another.  The  unknown  heat  transfer  coefficients 
that  are  being  determined  are  from  the  free  stream  to  the  fin 
area  of  the  vane,  not  the  mount  or  shaft.  The  small 
temperature  change  in  the  experimental  profiles  indicated  that 
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very  little  energy  from  the  fin  section  was  being  transmitted 
to  the  shaft  and  mount.  The  unknowns  were  identified  as  bn, 
b22/  b32,  b42  and  The  assun^Jtions,  for  the  first  series  of 
simulations,  were  that  b^  was  a  constant  and  b22=b32=b42.  With 
the  program  modification  in  place  the  simulation  was  conducted 
using  the  ZXSSQ  optimizer.  The  result  of  this  simulation  was 
the  same  as  in  the  six  node  model,  the  program  failed  to 
converge  using  this  optimizer.  Attenpts  were  made  at  using 
only  a  single  experimental  temperature  profile,  but  these 
simulations  also  failed  to  converge.  The  five  node  program 
was  modified  as  was  the  six  node  to  use  a  constrained 
optimizer  to  place  a  set  of  boundary  limits  on  the  unknowns 
with  the  limits  set  to  the  same  values  as  in  the  six  node 
simulations.  The  assumptions  for  the  five  node  model  unknowns 
remained  the  same  with  constant  and  b22=b32=b43.  This 
simulation  yielded  the  same  results  as  the  previous  modeling 
attempts  in  that  it  failed  to  converge.  A  more  systematic 
approach  was  taken  to  determine  if  the  simulations  could  be 
made  to  converge.  To  accomplish  this  the  simulations  would  be 
constructed  to  match  only  one  temperature  profile  at  a  time. 
The  assumptions  were  that  none  of  the  unknowns  were  dependent 
on  one  another  and  the  term  bj,  was  no  longer  considered  a 
constant  as  before.  The  unknown  variables  were  identified  as 
^11  ♦  ^«2  and  a^g. 

The  simulations  converged  and  succeeded  in  providing  an 
adequate  match  for  the  selected  single  node  time  ten^erature 
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profiles.  The  other  two  fin  area  node  profiles  were 
unsatisfactory  in  the  match  between  the  experimental  and 
calculated  profiles.  The  model  was  modified  to  match  two 
tenperature  profiles  from  the  fin  area  of  the  vane.  This 
series  of  simulations  converged  and  an  adequate  match  was 
achieved  for  two  of  the  three  tenperature  profiles,  but  again 
the  problem  was  in  the  third  profile.  The  program  matched  the 
two  temperature  profiles  that  it  was  set  up  to  match  but  the 
calculated  profile  for  the  third  node  was  substantially  off 
from  the  experimental  data.  Being  able  to  match  two  node 
temperature  profiles  the  program  was  modified  in  an  attempt  to 
match  all  three  nodal  temperature  profiles.  The  program  was 
modified  again  to  read  experimental  data  for  nodes  two,  three 
and  four,  and  to  match  all  three  time  temperature  profiles. 
This  simulation  was  conducted  using  the  assumptions  that  b^, 
baa*  ^32,  b4j,  and  Sso  were  independent,  the  recovery 
temperatures  were  Tni  =2670  and  T,a  =2570  Kelvin,  and  the  initial 
conditions  set  to  zero.  This  simulation  converged  and  yielded 
the  profiles  displayed  in  Figure  3.. 17.  The  simulation  was 
tested  for  its  robustness  by  starting  with  various  initial 
conditions.  The  program  converged  to  the  same  set  of  values 
with  the  same  set  of  time-temperature  profiles.  The  match 
between  the  calculated  and  experimental  profiles  of  node  two 
were  satisfactory  while  the  match  for  nodes  three  and  four 
could  only  be  considered  adequate.  The  deficiency  here  is  not 
with  the  modeling  program  but  with  the  experimental  time- 
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teirperature  profiles.  The  experimental  data  does  not  follow 
an  exponential  curve,  where  as  the  modeling  program  is  a  first 
order  linear  model  and  can  only  match  exactly  an  exponential 
profile.  The  modeling  program  matches  these  profiles  the 
closest  it  can  to  an  exponential  curve.  Another  point  to  note 
is  the  calculated  teirperature  profiles  for  all  three  nodes  are 
almost  a  straight  line  with  very  little  curve  in  the 
calculated  values  in  any  of  these  three  temperature  profiles. 
Since  the  exact  recovery  tenperatures  were  unknown  the  three 
node  temperature  matching  simulation  was  conducted  using  two 
other  sets  of  profiles.  Table  4  lists  the  recovery 
tenperatures  used  and  the  corresponding  convective  heat 
transfer  coefficients. 

TABLE  4 


Recovery 

Temperatures 

T.,/T„ 

Convective  Heat  Transfer  Coefficient 
W/m*»K 

ht  h„3  Ki  Ki 

2670/2570 

2970/2870 

3550/3450 

22027  1150  1057  594 
23562  878  933  527 
36819  405  762  434 

The  values  of  the  heat  transfer  coefficients  displayed  in 
Table  4  are  of  the  same  magnitude  as  those  obtained  by  Oriels 
in  the  quarter  scale  simulations.  The  values  for  the  heat 
transfer  coefficients  from  quarter  scale  modeling  were 
heB30122  W/m^*K  and  hy=210  W/m*«K  (Ref.  3:p.  12]. 
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The  calculated  profiles  in  the  convergent  five  node  model 
simulation  were  nearly  straight  lines.  This  effect  was  caused 
by  the  reduction  of  the  conductive  rate  coefficients,  the 
values.  These  values  are  determined  from  equations  (2)  and 
(4),  the  resistance  and  capacitance  equations.  Consider  the 
area  A  as  a  length  squared  and  the  volume  V  as  a  length  cubed. 
The  value  of  is  equal  to  (RC)*S  and  inserting  the 
equivalent  values  from  equations  (2)  and  (4),  and  considering 
the  physical  properties  as  constants,  the  value  of  RC  becomes 
a  function  of  the  length  squared.  This  means  that  the  value 
of  ajj  is  a  function  of  the  (length^) ■^.  As  the  quarter  scale 
model  is  increased  in  size  the  conductive  coefficients  aij 
decrease  as  the  reciprocal  of  the  length  squared,  if  the  same 
idea  is  applied  to  the  convective  coefficients  it  will  be  seen 
that  the  value  of  RC  changes  as  the  length  changes  and  the 
convective  rate  coefficient,  bij,  changes  as  the  (length) 

As  the  model  goes  to  full  scale,  the  geometric  properties 
become  larger  and  thus  the  rate  coefficients  loecome  smaller. 
Since  the  conductive  values  decrease  as  the  reciprocal  of  the 
square  of  a  length,  where  the  convective  values  only  decrease 
as  the  reciprocal  of  the  length,  the  convective  rate 
coefficients  become  dominant. 

The  quarter  scale  modeling  conducted  by  Oriels  (Ref.  3], 
demonstrated  that  the  value  of  the  convective  rate 
coefficient,  b^,  from  the  tip  section  of  the  vane  could  be 
adjusted  with  minimal  effect  on  the  other  unknown  variables. 
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This  was  due  to  the  ablative  cooling  of  the  tip  due  to  the 
erosion  of  the  tip  area  of  the  vane  by  the  rocket  exhaust.  In 
the  full  scale  experimental  rocket  firings  there  was  minimal 
vane  erosion  as  cortqpared  to  the  quarter  scale  vane  tests .  The 
energy  in  the  full  scale  tests  was  not  being  carried  away  by 
the  ablative  cooling  but  was  being  transmitted  into  the  vane. 
This  was  one  of  the  primary  reasons  that  led  to  making  the 
value  of  the  convective  rate  coefficient,  bn,  a  variable  in 
the  full  scale  model  instead  of  the  constant  that  was  assumed 
based  on  the  quarter  scale  model  results. 
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CONCLUSIONS 


•  The  four  node  quarter  scale  model  cannot  be  directly 
scaled  up  to  a  full  scale  model  due  to  the  problems 
encountered  with  substantial  juii¥>s  in  the  unknowns  asg  and 
ags,  corresponding  to  the  shaft  and  mount  locations. 

•  The  use  of  a  constrained  optimizer  is  necessary  to  hold 
the  value  of  the  unknowns  to  positive  values  and  prevent 
the  equation  solver  from  experiencing  an  arithmetic  fault 
or  causing  a  fatal  error  because  the  energy  balance 
equations  became  stiff. 

•  The  value  obtained  in  the  quarter  scale  vane  modeling 
process  more  adequately  reflects  the  average  convective 
heat  transfer  coefficient  of  the  fin  area  where  as  the 
v,^l*aes  obtained  from  the  full  scale  modeling  display  a 
3iore  accurate  representation  of  the  heat  transfer  process 
in  the  vane. 

•  The  calculated  profiles  of  Figure  3.17  indicate  that  as 
the  vane  size  increases  the  convective  heat  transfer 
properties  become  the  dominant  player  in  the  heat  transfer 
process  of  the  vane. 

•  The  rise  in  the  shaft  and  mount  experimental  temperature 
profiles  is  likely  due  to  the  rocket  exhaust  plume 
impinging  directly  on  the  shaft  of  the  vane. 

•  Hie  value  of  the  convective  rate  coefficient,  bj,,  must  be 
considered  an  unknovm  variable  in  the  time* temperature 
matching  process,  due  to  the  minimal  amount  of  vane 
erosion  and  the  lessening  effect  of  ablative  cooling  in 
the  full  scale  vane. 


73 


FCai  FORIHER  STUDY 


REC< 


•  It  is  recommended  that  further  study  modify  the  five  node 
model  by  reducing  it  to  three  nodes  and  coir^aring  the 
results  to  see  if  a  single  node  in  the  fin  area  of  the 
vane  could  produce  an  acceptable  convective  heat  transfer 
coefficient. 

•  The  full  scale  model  needs  to  be  able  to  take  into  account 
the  effects  of  radiation  from  the  casing  of  the  rocksst 
firing  chamber.  Itiis  may  condensate  for  the  non- 
linearities  in  the  experimental  temperature  profiles. 

•  Vane  erosion  was  evident  in  the  quarter  scale  testing  and 
some  vane  erosion  occurred  in  the  full  scale  testing.  The 
model  needs  to  be  upgraded  to  include  a  conponent  which 
can  compensate  for  the  erosion  based  on  the  metal  content 
of  the  rocket  exhaust. 

•  Since  the  convective  rate  coefficient  cannot  be  considered 
a  constant  at  the  tip,  a  thermocouple  attached  to  the  tip 
of  the  vane  to  collect  the  experimental  time-tenperature 
data  would  greatly  enliance  the  thermal  model  deveiepanent . 

•  The  process  for  the  attachm^t  of  theoswcouples  needs  to 
be  addressed.  The  experimental  data  had  glitches  or 
spikes  in  the  temperature  profiles  during  various  times  of 
the  rocket  firing  indicating  that  the  therawcouple 
temporarily  lost  contact  with  the  vana. 
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APPEHDIX  A.  PIO  PROGRAM 

The  progrcun  contained  is  the  FORUL^N  code  main  program 
used  in  the  vane  modeling  process.  Hie  programs  name  was 
changed  to  NODES  to  distinguish  it  from  the  PID  program  using 
the  six  node  modeling  approach 
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Ri  *r>  Ri  If  Jf> 


Program  NODES 

This  program  is  the  PID  program  for  the  five  node  vane  model, 
external  temp 

integer  m,n,iparm(6) ,ibtype,ldf jac 
parameter  (m=183,n=5,ldf jac=m) 
real*8  rparm(7) ,x(n) ,f (ra) ,xjac(m,n) ,xg(n) , 
xlb(n) ,xub(n) ,xscale(n) ,fscale(m) , 
ssq, float, 

a(5,5) ,b(5,5) ,u(5) ,t2(61) ,t3(61) ,t4(61) ,ys(5,61) , 
rho , k , op , vt , vf , vs , atf , af  s , Itf , If  s , sf , 
atf 1 , atf 2 , atf 3 , vf 1 , vf 2 , vf 3 , ubl , ub2 
intrinsic  float 

common/datal/a , b, u, t2 , t3 , t4 , ys 

Open  files  for  data  input/output 

open ( 10 , name= • result . dat ' ,  status= ' new ' ) 
open ( 9 , name= ' temp . mat ' ,  status= ' new ' ) 
open ( 8 , name= ' da tarn . dat ‘ ,  status®  ^  old ' ) 
open ( 7 , name® • input . dat ' ,  status® 'old') 

read  in  experimental  data 

do  i®l,61 

read(8,*)  t2(i) 

enddo 
do  i®l,61 

read(8,*)  t3(i) 

enddo 
do  i®l,61 

read (8,*)  t4(i) 

enddo 
close (8) 

read  in  input  data 


rsad(7, *) 
read (7,*) 
rcad(7, *) 
read (7,*) 
read(7,*) 
read ( 7 , * ) 
read(7,*) 
read ( 7 , * ) 
read (7,*) 
read (7,*) 
read (7,*) 
read (7,*) 
read (7,*) 
read (7,*) 
read(7, *) 
read(7,*) 
read(7, *) 
close (7) 


rho,k,cp 


vt,  vfl,  vf2,  vf3,  vs 


atfl,  atf 2,  atf 3,  afs 

Itf,  ifs 


sf,  ubl,  ub2 
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full  scale  data 


rl2=^100.  o*ltf/  (k*atf  1) 
r23=100.0*ltf/(k*atf2) 
r34=100. 0*ltf/ (k*atf3) 
r45=100. 0*lfs/ (k*afs) 
cl=rho*cp*vt*0 . 000001 
c2=rho*cp*vf 1*0 . 000001 
c3=rho*cp*vf 2*0 . 000001 
c4=rho*cp*vf3*0. 000001 
c5=rho*cp*vs*0 . 000001 


scaled  data 

rl2=rl2/sf 

r23=r23/sf 

r34=r34/sf 

r45=r45/sf 

cl«cl*sf**3 

c2=c2*sf**3 

C3=c3*sf**3 

C4=c4*sf **3 

c5=c5*sf**3 

rfl=0.5 

do 

u(i)»0.0 

do  j=l/5 
a(i,j)-o.o 
b(i, j)»0.0 
enddo 

enddo 

a(l,2)-l/(cl*rl2) 

a(2,l)-l/(c2*rl2) 

a(2,3)»l/(c2*r23) 

a(3,2)-l/(o3*r23) 

a(3,4)"l/(c3*r34) 

a(4,3)“l/(c4*r34) 

a(4,5)»l/(c4*r45) 

a(5,4)»l/(c5*r45) 

b(l,l)“l/(cl*rfl) 

wirite(6,*)^  al2  a21  a23  a32' 

writet6,3002)  a(l,2) ,a(2, 1) ,a(2,3) ,a(3,2) 
write(6,*)'  a34  a43  a45  a54' 

write(6,9003)  a(3,4) ,a(4,3) ,a(4,5) ,a(5,4) 

002  format(5f9.3) 

303  format(4f9.3) 

a5g  "O.o 
b(l,i)-o.o 

b(2|2)«0.0 
b(3,2)«0.0 
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a(l,l)=-(a(l,2)+b(l,l)) 
a(2,2)=-(a(2,l)+a(2,3)+b(2,2) ) 
a(3,3)=-(a(3,2)+a(3,4)+b(3,2)) 
a(4,4)=-(a(4,3)+a(4,5)+b(4,2)) 
a(5,5)=-(a(5,4)+a5g) 

u(l)=ubl 

u(2)=ub2 

xg(l)=a5g 

xg(2)=b(l,l) 

xg(3)=b(2,2) 

xg(4)=b(3,2) 

xg(5)=b(4,2) 

set  up  parameters  for  DBCLSF  call 
do  lc=l,n 
xscale(]c)=1.0 
xlb(k)«0.0001 
end  do 

xub(l) =100.0 
xub(2)=100.0 
xub(3) =100.0 
xub(4) =100.0 
xub(5)=100.0 

do  1=1, m 
fscale(l)»l. 0 
end  do 

ibtype»0 


call  dbclsf ( temp , m , n , xg , Ibtype , xlb , xub , xscale , f scale , 
&  iparm,rparm,x,f ,xjac,ldf jac) 

calculate  unknown  resistances  and  capacitances 

aSg  »x(l) 
b(l,l)-x{2) 
b(2,2)«x(3) 
b(3,2)-x(4) 
b(4,2)=x(5) 


rfl  •l/(b(l,l)*cl) 
rf22-l/(b(2,2)*c2) 
rf23-l/(b(3,2)*c3) 
rf24»l/(b(4,2)*^C4) 
r5g«l/ (a5g*c5) 

print  and  save  results 

write(6,*)  •  a5g  bll  b22  b32' 

write(6,9009)  x(l) ,x(2) ,x(3) ,x(4) ,x(5) 


09  format(5fl2.2) 
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;UU 


rormat;(5ti2 .4) 


write(l0,*) 'a5g  b(l,l)  b(2,2)  b(3,2)' 

write(10,9000)  x(l),x(2),x(3),x(4),x(e) 
write(10, *) 
write(10, *) 

write(10, *) 'rfl  rf22  rf23  rf24  r5g' 

write ( 10 , 9000)  rf Irf 22 , rf 23 , rf 24 , r5g 

write  the  temp-time  data  for  MATLAB  analysis 

do  i=l,61 
tt=0.0768*float(i) 

write(9,9001)tt,ys(2,i)  ,ys(3 , i)  ,ys(4, i) ,t2 (i) , t3 (i) , t4 (i) 
enddo 

01  format(lf6.2,6fl0.3) 

close(lO) 

close(9) 

end 


Subroutine  TEMP  (m,n,x,f) 

This  calculates  the  temperature-time  history  using  the  current 
parameters  supplied  by  ZXSSQ  called  from  PID.  It  calculates 
an  error  function  returned  to  ZXSSQ  based  on  the  differences 
between  predicted  and  observed  temperature  histories. 

integer  maxparam,  neq 
parameter (maxparam^SO,  neq=5) 
integer  idO# istep,nout,m,n 

real*8  t,y(neq) ,a(neq,neq) ,b(neq,neq) ,u(neq) ,param(raaxparam) 
real*8  tend,tol, fen, float,a5g,ys(5, 61) ,x(n) , 

&  f (m) ,t2(61) ,t3(61) ,t4(61) 

intrinsic  float 
external  fcn,divprk,sset 

common/datal/a , b, u , t2 , t3 , t4 , ys 

open  ( 12 ,  name®  •  incoming .  dat  • ,  .status® '  new  • ) 

aSg  »x(l) 
b(l,l)»x(2) 
b(2,2)-x(3) 
b(3,2)«x(4) 
b(4,2)»x(5) 

write(6,8000)a5g,b(l,l) ,b(2,2),b(3,2),b(4,2) 

00  format(5fl2.4) 

a(l,l)  — (a(l,2)+b(l,l)) 
a\'2,2)  — (a(2,l)+a(2,3)+b(2,2) ) 
a(3,3)  — (a(3,2)+a(3,4)+b(3,2)) 
a(4,4)  — (a(4,3)+a(4,5)+b(4,2)) 
a(5,5)  — (a(5,4)+a5g) 
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t=0.0 

do  i=l,neq 
y(i)=»0.0 

do  j=l,6l 
ys(i, j)=0.0 
enddo 

enddo 

tol=0.0005 

call  sset  (maxparm,  0.0«  param,  1} 


idO=l 

do  isteps'ljei 

tend=0.0768*float (istep) 

CALL  DIVPRK  (idO,  neq,  fen,  t,  tend,  tol,  param,  y) 

do  i=l,5 

ys(i,istep)=y(i) 

enddo 

enddo 

Final  call  to  release  workspace 
id0=3 

call  divprk  (idO, neq, fen, t, tend, tol, parain,y) 

calculate  error  functions 

do  i=l,61 
f(i)-ys(2,i)-t2(i) 
f  (i+6l)-y8(3,i)-t3(i.) 
f (i+122)»y8(4,i)-t4(i) 
enddo 

print  out  rms  error 
S8qr»o . 0 
do  i-1,122 
ssqr»ssqr +f ( i ) *  f ( i ) 
enddo 

ssqr»ssqr/122 

xer“sqrt(ssqr) 

write (6,*)  xer 

write (12,*)  xer 

return 

end 


subroutine  fcn(neq,t,y,yprime) 
integer  neq 

real*8  t,y(neq)  ,ypriine(neg) 
real*8  a(5,5) ,b(5,5) ,u(5) ,d,ys(5,61) 
conmon/datal/a , b, u , t2 , t3 , t4 , ys 

thrust  profile  simulation  as  step  input 


if  (t.gt.0.2)  then 
d*1.0 
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else 

d=0.0 
end  if 

do  i=l,neq 
ypriine(i)  =0.0 

do  j=l,neq 

ypriine(i)=yprime(i)+a(i,  j)*y(j)+b(i,  j)*u(j)  *d 
enddo 

enddo 

return 

end 
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APPENDIX  B.  GEOMETRIC  AND  PHYSICAL  DATA  FILE 


The  file  in  this  appendix,  INPUT.DAT,  is  the 
containing  the  geometric  and  physical  properties  of  the 
The  values  are  listed  in  full  scale. 


file 

vane. 
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NAWC  INVERSE  HEAT  TRANSFER  PROGRAM.  INPUT  DATA  FOR  FULL  SCALE 

Material  properties:  rho  (kg/in''3) ,  k  (w/m-deg  k) ,  Cp  (J/kg  deg  k) 

18310.0,  173.0,  146.0 

Vol  (tip,  cin*3) ,  Vol  (fin,  cm*3) ,  Vol  (shaft,  cm^3) 

2.6  12.86  17.15  21.44  23.0 

X-section  areas:  tip-fin  (cm'' 2)  fin-shaft  (cm'' 2) 

5.6  5.9  6.2  5.2 

Conductive  path  tip-fin  (cm)  fin-shaft  (cm) 

2.774  5.0 

Scale  factor: 

1.00  2670  2570 
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APPENDIX  C.  TPROFILE  -  FORNARD  MODELING  PROGRAM 
The  program  this  appendix  is  the  forward  modeling  program 
based  on  a  modification  of  the  main  program,  PID,  and  used  to 
illustrate  the  effects  of  the  unknown  variables  on  the 
resultant  time  tenperature  profiles. 
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ta 


Program  TPROFILE 


This  is  a  forward  model  program  developed  from  the  program 
PID  by  removing  the  optimizer. 

external  temp 

integer  m,n, ixjac,nsig,maxfn, iopt, infer, ier 
parameter  (m=183 , n^S) 

real*8  parm(6) ,x(n) ,f (m) ,xjac(m,n) /Xjtj { (n+1) *n/2) , 
work(5*n+2*m+( (n+l) *n/2) ) , 
eps , de 1 ta , ssg , ubl , ub2 , 

a(5,5) ,b(5,5) ,U(5) ,t2(61) ,t3(61) ,t4(61) ,ys(5,6?)/ 
rho,k,cp,vt, vf ,vs,atf ,afs, Itf , lfs,sf 
common/datal/a , b, u , t2 , t3 , t4 , ys 

Open  files  for  data  input/output 

open (10, name- 'result. dat ' ,  s t atus= ' new ' ) 
open ( 9 , name® ' temp . mat ' ,  status* ' new ' ) 
open ( 8 , name* ' datam . dat ' ,  status* ' old ' ) 
open ( 7 , name* ' input . dat ' ,  status* • old ' ) 

do  i*l,6l 

read (8,*)  t2(i) 

enddo 
do  i“l,61 

read(8,*)  t3(i) 

enddo 
do  i*l,61 

read(8,*)  t4(i) 

enddo 
close (8) 

read  in  input  data 

read (7,*) 
read(7,*) 
read (7,*) 
read(7, *) 

read (7,*)  rho,k,op 

read(7,*) 

read(7,*) 

read (7,*)  vt,  vfl,  vf2,  vf3,  vs 

read(7,*) 

read  (7, 

read (7,*)  atfl,  atf2,  atf3,  afs 
read (7,*) 
read (7,*) 

read(7,*)  Itf,  Ifs 

read(7,*) 

read(7,*) 

read (7,*)  sf,  ubl,  ub2 
cloae(7) 

initial  conditions 


Cull  scale  data 
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rl2=l00 . 0*ltf / (k*atf 1) 
r23=100.0*ltf/ (k*atf2) 
r34=l00 . o*ltf / (k*atf 3) 
r45=100.0*lfs/ (k*afs) 
cl=rho*cp*vt*0 . 000001 
c2=rho*cp*vf 1*0 . 000001 
c3=rho*cp*vf 2*0. 000001 
c4=rho*cp*vf 3*0. 000001 
c5=rho*cp*vs*0 . 000001 

write ( 6  ^  *) rl2 , r23 , r34 , r45 , cl , c2  ^  c3 , c4 , c5 

scaled  data 

rl2=rl2/sf 

r23=r23/sf 

r34=r34/sf 

r45=r45/sf 

cl=cl*sf**3 

c2=c2*sf**3 

c3=c3*sf**3 

c4®c4*sf**3 

C5»c5*sf**3 

rf 1-0. 076318 

do  i-1,5 
u(i)»0.0 

do  j-l,5 
a(i, j)=0.0 
b(i, j)«0.0 
enddo 

anddo 

a(l,2)-l/(cl*rl2) 

a(2,l)-l/(c2*rl2) 

a{2,3)-l/(c2*r23) 

a(3,2)»l/(c3*r23) 

a(3,4)-l/(c3*r34) 

a(4,3)-l/(C4*r34) 

a(4,5)-l/(c4*r45) 

a(5,4)«l/(c5*r45) 

b(l,l)»j./(cl*rfl) 

write(6,*)'  al2  a2l  a23  a32  bll' 

write(6,9002)  a(l,2) ,a(2,l) ,a(2,3) ,a(3,2) ,b(l,l) 
write(6,*)'  a34  a43  a45  a54' 

write(6,9003)  a(3,4) ,a(4,3) ,a(4,5) ,a(5,4) 

)02  format(5f9.3) 

)03  Corinat(4f9.3) 

aSg  -0.0001 

b(l,l)-300.0 

b(2,2)-1.339 

b(3,2)-0.56S7 

b(4,2)-0.0037 
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a\j.,  J./ - vcna.,6/-rw^x,x;; 

a(2,2)=-(a(2,l)+a(2,3)+b(2,2)  ) 
a(3,3)=-(a(3,2)+a(3,4)+b(3,2)) 
a(4,4)=-(a(4,3)+a(4,5)+b(4,2)) 
a(5,5)=-(a(5,4)+a5g) 

u(l)=ubl 

u(2)=ub2 

x(l)=a5g 

x(2)»b(l,l) 

x(3)=b(2,2) 

x(4)=b(3,2) 

x(5)=b(4,2) 

call  teinp(x,m,n, f ) 


wrice  the  temp>time  data  for  MATLAB  analysis 
do 

tt».0768*float(i) 

write{9,900l)tt,ys(2,i) ,ys(3,i) ,ys(4,i) ,t2(i) ,t3(i) ,t4(i) 
enddo 

01  £orinat(lf6.2,6fl0.3) 

close (10) 
close (9) 
end 


Subroutine  TEMP 

This  calculates  the  temperature-time  history  using  the  current 
parameters  supplied  by  ZXSSQ  called  from  PID.  It  calculates 
an  error  function  returned  to  ZXSSQ  based  on  the  differences 
between  predicted  and  observed  temperature  histories. 

integer  maxparam,  neq 
parameter (maxparam^SO,  neq«S) 
integer  idO,istep«nout^m,n 

real*8  t,y(neq) ,a(neq,neq} «b(neq,neq) ,u(neq) ,param(maxparam) 
real*8  tend, tol, fen, float,a5g,ys( 5,61} ,x(n) , 

&  f (m) ,t2(61) ,t3(61) ,t4(61) 

intrinsic  float 
external  fen , di vprk , sset 

common/ da ta l/a,b,u,t2,t3,t4,ys 

aSg  -x(l) 
b(l,l)-x(2) 
b(2.2}»x(3) 
b(3,2)»X(4) 
b(4,2)-X(5) 

a(l,l)*-(a(l,2)+b(l,l)) 
a(2,2)  — (a(2,l)+a(2,3)+b(2,2)) 
a(3,3)  — (a(3,2)>a(3,4)+b(3,2)) 
a(4,4)  — (a(4,3)+a(4,5)-<-b(4,2)) 
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Set  initial  conditions 


t=0.0 

do  i=l,neq 
y(i)=0.0 

do  j=l,61 
ys(i, j)=0.0 
enddo 

enddo 

tol=0.0005 

call  sset  (maxparm,  0.0,  par am,  1} 
param(4)=2000 


id0=l 

do  istep=l,61 

tend=0 . 07  68  loat { istep) 

CALL  DIVPRK  (ido,  neq,  fen,  t,  tend,  tol,  param,  y) 

do  is»l,5 

ys(i, istep) «y(i) 

enddo 

enddo 

Final  call  to  release  workspace 


id0a3 

call  divprk  ( idO, neq, fen, t, tend, tol, param,y) 

do  i-l,61 

f(i)-ys(2,i)-t2(i) 

f(i)-ys(3,i)-*t3(i) 

f(i)-ys(4,i)-t4(i) 

enddo 

ssqr-O.O 
do  1^1,61 

ssqr»&sqr+f ( i ) *  f ( i ) 
enddo 

ssqr»ssqr/6l 

xer»sqrt(S8qr) 

write(6,*)  xer 

return 

end 


subroutine  fcn(neq,t,y,yprinie} 
integer  neq 

realms  t,y(neq) ^yprime(neq) 

real*8  a(5, 5) ,b(5, 5) ,u(5) ,d,ys(5,6l) 

CQmmon/datal/a,b,u,t2,t3,t4,ys 

if  (t.gt.0.2)  then 
d-l.O 

else 


d»0.0 


88 


do  i=l,neq 
yprime(i)=0.0 

do  j=l,neq 

ypriine(i)=yprinie(i)+a(i,  j)*y(j)+b(i,  j)*u(j)*d 
enddo 

enddo 


return 

end 


APPEKIDIX  D.  CONSTRAINED  OPTIMIZER  -  DBCLSF 


This  optimizer  is  the  constrained  optimizer  which  allows 
the  user  to  place  constraints  on  the  range  and  values  of  the 
unknowns.  The  constraints  may  be  a  set  of  boundary  limits  or 
restraining  the  unknown  variables  to  either  being  negative  or 
positive. 
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BCLSJ/DBCLSJ  (Single/Double  precision) 

Purpose:  Solve  a  nonlinear  least  squares  problem  subject  to  bounds 

on  the  variables  using  a  modified  Levenberg-Marquardt 
algorithm  and  a  user-supplied  Jacobian. 

Usage:  CALL  BCLSJ  (FCN,  JAC,  M,  M.  XGUESS,  IBTYPE,  XLB.  XUB, 

XSCALE,  FSCALE,  IPARAM,  RPARAM,  X.  FVEC, 

FJAC,  LDFJAC) 

Arguments 

FCN  -  User-supplied  SUBROUTINE  to  eveduate  the  function  to  be 
minimized.  The  usage  is 
CAa  FCH  (M.  M.  X.  F).  where 

M  -  Length  of  F.  (Input) 

N  -  Length  of  X.  (Input) 

X  -  The  point  at  which  the  function  is  evaluated. 
(Input) 

X  should  not  bo  changed  by  FCN. 

F  «  The  computed  function  at  the  point  X. 

(Output) 

FC?I  must  be  declared  EXTERNAL  in  the  calling  program. 

JAC  -  User-supplied  SUBROUTINE  to  evaluate  the  Jacobian  at  a 
point  X.  The  usage  is 
CALL  JAC  (M.  N,  X.  FJAC,  LDFJAC).  where 
M  -  Length  of  F.  (Input) 

N  -  Length  of  X.  (Input) 

X  -  The  point  at  which  the  function  is  evaluated. 
(Input) 

X  should  not  be  changed  by  FCN. 

FJAC  -  The  computed  M  by  N  Jacobian  at  the  point  X. 
(Output) 

LOFJAC  -  Leading  dimension  of  FJAC.  (Input) 

JAC  must  be  declared  EXTERNAL  in  the  calling  program. 

N  -  Number  of  functions.  (Input) 

N  -  Number  of  variables.  (Input) 

XGUESS  -  Vector  of  length  N  containing  the  initial  guess.  (Input) 
IBTYPE  -  Scalar  Indicating  the  types  of  bounds  on  variables. 

(Input) 

IBTYPE  Action 

0  User  will  supply  all  the  bounds. 

1  All  variables  are  nonnegatlve. 

2  All  variables  are  nonpositive. 

3  User  supplies  only  the  bounds  on  1st  variable, 
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all  other  variables  will  have  the  same  bounds. 

XLB  -  Vector  of  length  N  containing  the  lower  bounds  on 

variables.  (Input,  if  IBTYPE  *  0;  output,  if  IBTYPE  »  1 
or  2;  input/output,  if  IBTYPE  »  3) 

XUB  Vector  of  length  M  containing  the  upper  bounds  on 

variables.  (Input,  if  IBTYPE  ■  0;  output,  if  IBTYPE  »  1 
or  2;  input/output,  if  IBTYPE  *  3) 

XSCALE  -  Vector  of  length  N  containing  the  diagonal  scaling  matrix 
for  the  variables.  (Input) 

In  the  absence  of  other  information,  set  all  entries  to 

1.0. 

FSCALE  -  Vector  of  length  M  containing  the  diagonal  scaling  matrix 
for  the  functions.  (Input) 

In  the  absence  of  other  information,  set  all  entries  to 

1.0. 

IPARAM  -  Parameter  vector  of  length  6.  (Input/Qutput) 

See  Remarks. 

RPARAM  ~  Parameters  vector  of  length  7.  (Input/Output) 

See  Remarks. 

X  -  Vector  of  length  N  containing  the  approximate  solution. 
(Output) 

FVEC  -  Vector  of  length  H  containing  the  residuals  at  she 
approximate  solution.  (Output) 

FJAC  •  H  by  M  matrix  containing  a  finite  difference  approximate 
Jacobian  at  the  approximate  solution.  (Output) 

LOFJAC. •  Leading  dimension  of  FJAC  exactly  as  specified  in  the 
dimension  stattiaent  of  the  calling  program.  (Input) 

Remarks 

L.  Automatic  workspace  usage  is 

BCLSJ  14«N  2*H  *  I  units,  or 

08CLSJ  26*N  ♦  4*M  -  2  units. 

Uorkspaee  may  be  explicitly  provided,  if  desired,  by  use  of 
B2LSJ/0B2LSJ .  The  reference  is 

CAa  B2I.SJ  (FCN,  JAC.  M,  H.  XCUESS,  IBTYPE,  XLB.  XUB. 

XSCALE.  FSCALE.  IPARAM,  RPARAN.  X.  FVEC. 

FJAC.  LOFJAC.  WK.  IVK) 

The  additional  arguments  are  as  follows: 

WX  -  Work  vector  of  length  12«N  *-  2«H  -  1.  WK  contains 
the  following  information  on  output; 

The  second  K  locations  contain  the  last  step  taken. 

The  third  locations  contain  the  last  Gauss^Newton  step. 
The  fourth  M  locations  contain  an  estimate  of  the 
gradient  at  the  solution. 
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IWK  -  Work  vector  of  length.  2«N  containing  the 

permutations  used  in  the  QR  factorization  of  the  Jacobian 
at  the  solution. 

2.  Informational  errors 
Type  Code 

3  1  Both  the  actual  and  predicted  relative  reductions  in  the 

function  are  less  than  or  equal  to  the  relative  function 
convergence  tolerance. 

4  2  The  iterates  appear  to  be  converging  to  a  noncritical 

point. 

4  3  Maximum  number  of  iterations  exceeded. 

4  4  Maximum  number  of  function  evaluations  exceeded. 

4  5  Five  consecutive  steps  have  been  taken  with  the  maximum 

step  length. 

4  6  Maximum  number  of  Jacobian  evaluations  exceeded. 

3.  The  first  stopping  criterion  for  BCLSJ  occurs  when  the  norm  of  the 
function  is  less  than  Che  absolute  function  colerancs.  The  second 
stopping  criterion  occurs  when  the  norm  of  the  scaled  gradient  is 
less  chan  the  given  gradient  tolerance.  The  third  stopping 
criterion  for  BCLSJ  occurs  when  the  scaled  distance  between  the 
last  two  steps  is  less  than  the  step  tolerance. 

4.  If  nondcfault  parameters  are  desired  for  IPARAm  or  RPARAK,  then 
U4LSF  is  called  and  the  corresponding  parameters  are  set  to  the 
desired  value  before  calling  the  optimization  program.  Otherwise, 
if  Che  default  parameters  are  desired,  then  set  IPARAH(l)  to  zero 
and  call  the  optimization  program  omitting  the  call  to  U4LSF. 

The  call  to  OdLSF  would  be  as  follows: 

CALL  U4LSF  (IPARAN,  aPARAM). 

The  following  is  a  list  of  the  parameters  and  the  default  values: 
IPARAM  •  Integer  vector  of  length  6. 

IPARAM(l)  Initialization  flag.  (0) 

IPARAN(2)  •  Mumber  of  good  digits  in  the  function. 

(Machine  dependent) 

IPARANO)  ■  .Maximum  number  of  iterations.  (100) 

IPARAN(4}  ■  Maximum  number  of  function  evaluations.  (400) 
IPARAM(S)  ■  Maximum  number  of  Jacobian  evaluations.  (100) 
IPARAN(6)  «  Internal  variable  scaling  flag.  (1) 

If  IPARAM(6)  >  I  Che  values  for  ISCALE  are 
set  internally. 

aPARAM  ■*  Real  vector  of  length  ?. 
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RPARAM(l)  «  Scaled  gradient  tolerance.  (eps*«(2/3)) 
RPAilAM(2)  •  Scaled  step  tolerance.  (eps*»(2/3)) 

RPARAM(3)  ■  Relative  function  tolerance. 

(MAX ( 1 . OE-10 . eps»* (2/3) ) ) 

RPARAM<4)  «  Absolute  function  tolerance. 

(MAX( 1 . OE-10 , eps** (2/3) ) ) 

RPARAM(S)  »  False  convergence  tolerance.  (100*eps) 
RPARAM(6)  •  Maximum  allovable  step  size. 

(1000-MAX (ton. T0L2))  where. 

TQLl  •  SQRKsum  of  (XSCAlEO-XGUESSd))— 2) 
for  I  ■  l, . . . ,M 
TQL2  •  2-norm  of  XSCALE. 

RPARAM(7)  «  Size  of  initial  trust  region  radius. 

(Based  on  the  initial  scaled  Cauchy  stop) 
eps  is  machine  epsilon. 

If  double  precision  is  desired,  then  0U4LSF  is  called  and  RPARAM 
is  declared  double  precision. 

Keywords:  Levenberg-Marquardt ;  Trust  region 


Algorithm 

3CLSJ  uses  a  modified  Levenberg<Marquardt  method  and  an  active  set  strategy  to 
solve  nonlinear  least  squares  problems  subject  to  simple  bounds  on  the  variables. 
The  problem  is  stated  as  follows: 

mmjFU)^F{r)=5f;A(x)= 

subject  to  /  <  X  <  u. 

where  m  >  n.  F  :  R**  —  R"*.  and  /,{x)  is  the  t-th  component  function  of  F(x). 
From  a  given  starting  point,  an  active  set  lA.  which  contains  the  indices  of  the 
variables  at  their  bounds,  is  built.  variable  is  called  a  free  variable'  if  it  is  not  in 
the  active  set.  The  routine  then  computes  the  search  direction  for  the  foee  variables 
according  to  the  formula 

da 

where  p  is  the  Leveoberg*Marquardt  parameter.  F  »  F(x).  and  J  is  the  Jacobian 
with  respect  to  the  free  variables.  The  search  direction  for  the  variables  in  lA  is 
set  to  zero.  Ttie  trust  region  approach  discussed  by  Dennis  and  Schnabel  (19d3) 
is  used  to  find  the  new  point.  Finalty.  the  optimaiity  conditions  iue  checked.  The 
conditions  are: 

Ifoi-c.HI  <  <•  /.  <x,  <f*, 
yix.i  <0.  r,  =  Ut 
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gix^)  >0.  ar,  =  li, 

where  e  is  a  gradient  tolerance.  This  process  is  repeated  until  the  optimality  criterion 
is  achieved. 

The  active  set  is  changed  only  when  a  free  variable  hits  its  bounds  during  an 
iteration,  or  the  optimality  condition  is  met  for  the  free  variables  but  not  for  adl 
variables  in  lA,  the  active  set.  In  the  latter  case,  a  variable  which  violates  the 
optimality  condition  will  be  dropped  out  of  lA.  For  more  detail  on  the  Levenberg- 
Marquardt  method,  see  Levenberg  (1944).  or  Marquardt  (1963).  For  more  detailed 
information  on  active  set  strategy,  see  Gill  and  .Murray  (1976). 


Example 

The  nonlinear  least  squares  problem 


min 


“  »al 


subject  to  -2  <  Xi  <  0.5 

-1  <-C2  <  2. 

where  /i(x)  »  I0(x2  -  xf).  and  /i(x)  =*  (I  -  xi)  is  solved  with  an  initial  guess 
(-1.2. 1.0).  and  default  values  for  parameters. 


C 

C 


C 

C 


C 

c 

c 


Oeclaracloa  of  variables 

ISTEQER  LOFJAC.  H.  H 
PAAAKETER  (LI)PJAC«2.  K-2.  .M-2) 


iKTEGEa  mKwa).  itp,  sour 

REAt  FJAC(10FJAC..1) .  FSCALEIH).  PVCC(M).  ROSBCX.  iUlSJAC. 

k  aPARAN(7).  X(M}.  XQUESSd).  XLa(8} .  IS(it).  XUB(.'I) 

eXTEUIAL  SCLSJ.  ROSSa.  ROSJAC.  tlNACH 

Coapate  the  least  squares  for  the 
RoseabFoch  fuact;oa 


DATA  XGUESS/-1.2E0.  l.OEO/.  XS/2-I.OEO/,  FSC' 

DATA  XL3/-2.0E0.  -l.OEO/,  XOB/O.SEO.  2.0E0/ 

All  the  bounds  are  (>rovided 


ITP  -  0 


IPARANd)  •  0 


Default  parsMters  are  used 


CAU.  BCLSJ  (R0SBCS.  ROSJAC.  M.  H.  {GUESS,  ITP,  XL3.  XUB.  IS, 
A  FSCALE,  IPARAM.  RPAAAN.  X.  FVEC.  FJAC.  LOFJAC) 

Priac  results 

CALL  UNACU  (2.  SODT) 

VRITE  ( ROUT. 09999)  X.  FVEC.  IPARANO).  IPARAM(4) 


99999  -ORNAT  '  The  solutloo  is  ’.  2F9.4.  //.  '  The  functj.oo  '. 
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i  ‘evaluated  at  the  solutioa  is  18X,  2F9.4,  //, 

4  '  The  auober  of  iteratioas  is  lOX,  13,  /,  '  The 

t  ‘niuBber  of  function  evaluations  is  13,  /) 

END 

C 

SUBROOTINE  ROSBCX  (M,  K,  X.  F) 

INTEGER  N,  N 

REAL  X(N),  F<M) 

C 

F(l)  -  1.0E1-(X(2)-X(1)-X(1)) 

F(2)  »  l.OEO  -  X(l) 

RETURN 

END 

C 

SUBROUTINE  RQSJAC  (M.  N.  X,  FJAC,  LOFJAC) 

INTEGER  M,  N.  LOFJAC 

REAL  X(N).  FJAC(LDFJAC.N) 

C 

FJAC(l.l)  •  -20.0E0*X(l) 

FJAC(2.1)  •  -l.OEO 
FJAC(l,2)  •  lO.OEO 
FJAC(2.2)  •  O.OEO 
RETURN 
END 


Output 


The  solutioa  is  .5000  .2500 

Th«  foaetioa  svaluatsd  at  th«  solution  is 
.0000  .5000 

Th«  auober  of  iteratioas  is  13 

the  number  of  function  evaluations  is  2i 
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APPENDIX  E.  ZXSSQ  OPTIMIZER 


This  optimizer  routine  was  the  initial  optimizing  program 
in  the  vane  modeling  process  using  the  IMSL  subroutines.  It 
provides  no  way  of  constraining  the  values  of  the  unknown 
variables . 
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I.MSL  ROUTINE 
PURPOSE 

USAGE 

ARCUHENTS 


NAME  -  ZXSSQ  ' 

<  MINIMUM  OP  the'  sum  OP  SQUARES  OP  M  PUNCTIONS 
IN  N  VARIABLES  USING  A  PINITE  OIFPERENCE 
LEVENBERG-HARQUAROT  ALGORITHM 

-  CALL  ZXSSQtPUNC, M,N,NSIC, EPS, DELTA, MAXPN, lOPT, 
P ARM, X.SSQ.F.XJAC.IXJAC.XJTJ, WORK, INFER. lER) 

PUNC  -  A  USER  SUPPLIED  SUBROUTINE  WHICH  CALCULATES 
THE  RESIDUAL. VECTOR  F( 1 ), F(2 F( M )  FOR 
GIVEN  PARAMETER  VALUES  X( 1 ). X( 2 X(  N) . 
THE  CALLING  SEQUENCE  HAS  THE  FOLLOWING  FORM 
CALL  PUNCCX.M.N.F) 

WHERE  X  IS  A  VECTOR  OF  LENGTH  N  AND  F  IS 
A  VECTOR  OP  LENGTH  M. 

FUNC  MUST  APPEAR  IN  AN  EXTERNAL  STATEMENT 
IN  THE  CALLING  PROGRAM. 

FUNC  MUST  NOT  ALTER  THE  VALUES  OF 
X(I).1-1 . N,  M,  OR  N. 

H  ^  THE  NUMBER  OF  RESIDUALS  OR  OBSERVATIONS 

( INPUT) 

N  -  THE  NUMBER  OF  UNKNOWN  PARAMETERS  (INPUT). 

NSIQ  -  FIRST  CONVERGENCE  CRITERION.  (INPUT) 

CONVERGENCE  CONDITION  SATISFIED  IF  ON  TWO 
SUCCESSIVE  ITERATIONS,  THE  PARAMETER 
ESTIMATES  AGREE,  COMPONENT  BY  COMPONENT. 

TO  NSIG  DIGITS. 

EPS  •  SECOND  CONVERGENCE  CRITERION.  (INPUT) 

CONVERGENCE  CONDITION  SATISFIED  IF.  ON  TWO 
SUCCESSIVE  ITERATIONS  THE  RESIDUAL  SUM 
OF  SQUARES  ESTIMATES  HAVE  RELATIVE 
DIFFCRENCE  LL.^iS  THAN  OR  EQUAL  TO  EPS.  CPS 
MAY  BE  SET  TO  YERO. 

DELTA  -  THIRD  CONVERGENCE  CRITERION.  (INPUT) 

CONVERGENCE  CONDITION  SATISFIED  IF  THE 
(EUCLIDEAN)  NORM  OF  THE  APPROXIMATE 
GRADIENT  IS  LESS  THAN  OR  EQUAL  TO  DELTA. 
DELTA  HAY  BE  SET  TO  ZERO. 

NOTE.  THE  ITERATION  IS  TERMINATED,  AND 
CONVERGENCE  IS  CONSIDCREO  ACHIEVCD.  IF 
ANY  ONE  OF  THE  THREE  CONDITIONS  IS 
SATISFIED. 

NAXFN  •  INPUT  MAXIMUM  NUMBER  OF  FUNCTION  EVALUATIONS 
(I.C..  CALLS  TO  SUBROUTINE  FUNC)  ALLOWED. 

THE  ACTUAL  NUMBER  OF  CALLS  TO  FUNC  MAY 
EXCEED  MAXFN  SLIGHTLY. 

lOPT  INPUT  OPTIONS  PARAMETER. 

lOPT-O  IMPLIES  BROWN'S  ALGORITHM  WITHOUT 
STRICT  DESCENT  IS  DESIRED. 
tOPT>l  IMPLIES  STRICT  DESCENT  AND  DEFAULT 
VALUES  FOR  INPUT  VECTOR  FARM  ARE  DESIRED. 
lOPt-2  IMPLIES  STRICT  DESCENT  IS  DESIRED  WITH 
USER  PARAMETER  CHOICES  IN  INPUT  VECTOR  FARM. 

PABM  '  INPUT  VECTOR  OF  LENGTH  N  USED  ONLY  FOR 
lOPT  EQUAL  TWO.  PARMd)  CONTAINS,  WHEN 

ZEROS  AND  EXTREMA  ZXSSQ- 1 
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!•) ,  THE  INITIAL  VALUE  OF  THE  MARQUAROT 
PARAMETER  USED  TO  SCALE  THE  DIAGONAL  OF 
THE  APPROXIMATE  HESSIAN  MATRIX.  XJTJ , 

BY  THE  FACTOR  (1.0  ♦  PARM(l)).  A  SMALL 
VALUE  GIVES  A  NEWTON  STEP,  WHILE  A  LARGE 
VALUE  GIVES  A  STEEPEST  DESCENT  STEP. 

THE  DEFAULT  VALUE  FOR  PARM(l)  IS  0.01. 

1-2.  THE  SCALING  FACTOR  USED  TO  MODIFY  THE 
MARQUAROT  PARAMETER,  WHICH  IS  DECREASED 
BY  PARH(2)  AFTER  AN  IMMEDIATELY  SUCCESSFUL 
DESCENT  DIRECTION.  AND  INCREASED  BY  THE 
SQUARE  OF  PARH(2)  IF  NOT.  PARM(2}  MUST 
BE  GREATER  THAN  ONE,  AND  TWO  IS  DEFAULT. 
1-3.  AN  UPPER  BOUND  FOR  INCREASING  THE 
MARQUAROT  PARAMETER.  THE  SEARCH  FOR  A 
DESCENT  POINT  IS  ABANDONED  IF  PARH(3)  IS 
EXCEEDED.  PARH(3)  GREATER  THAN  100.0  IS 
RECOMMENDED.  DEFAULT  IS  120.0. 

I-H.  VALUE  FOR  INDICATING  WHEN  CENTRAL 

RATHER  THAN  FORWARD  DIFFERENCING  IS  TO  BE 
USED  FOR  CALCULATING  THE  .HCOPIAN.  THE 
SWITCH  IS  MADE  WHEN  THE  NORM  OF  THE 
GRADIENT  OF  THE  SUM  OF  SQUARES  FUNCTION 
BECOMES  SMALLER  THAN  PARM(V),  CENTRAL 
DIFFERENCING  IS  GOOD  IN  THE  VICINITY 
OF  THE  SOLUTION.  SO  PARH(V)  SHOULD  BE 
SMALL.  THE  DEFAULT  VALUE  IS  0.10. 

X  -  VECTOR  OF  LENGTH  N  CONTAINING  PARAMETER 

VALUES. 

ON  INPUT.  X  SHOULD  CONTAIN  THE  INITIAL 
ESTIMATE  OP  THE  LOCATION  OF  THE  MlNIHUH. 

ON  OUTPUT.  X  CONTAINS  THE  FINAL  ESTIMATE 
OP  THE  LOCATION  OF  THE  MINIMUM. 

SSQ  -  OUTPUT  SCALAR  WHICH  IS  SET  TO  THE  RESIDUAL 
SUMS  OF  SQUARES.  F1 1 ) • *2 ♦ . . . • F(M ) * *2 ,  FOR 
THE  FINAL  PARAMETER  ESTIMATES. 

F  -  OUTPUT  VECTOR  OF  LENGTH  H  CONTAINING  THE 

RESIDUALS  FOR  THE  PINAL  PARAMETER  ESTIMATES. 

XJAC  -  OUTPUT  M  BY  N  MATRIX  CONTAINING  THE 

APPROXIMATE  JACOBIAN  AT  THE  OUTPUT  VECTOR  X. 

IXJAC  -  INPUT  ROW  DIMENSION  OF  MATRIX  XJAC  EXACTLY 
AS  SPECIFIED  IN  THE  DIMENSION  STATEMENT 
IN  THE  CALLING  PROGRAM. 

XJTJ  -  OUTPUT  VECTOR  OF  LENGTH  (N*l)*N/2  CONTAINING 

THE  N  BY  N  MATRIX  (XJAC-TRANSPOSCD)  •  (XJAC) 
IN  SYMMETRIC  STORAGE  MODE. 

WORK  -  WORK  VECTOR  OF  LENGTH  5«N  ♦  2*M  ♦  {N*l)»N/2. 

ON  OUTPUT*  WORK(I)  CONTAINS  FOR 

1-1.  THE  NORM  OP  THE  GRADIENT  DESCRIBED 
UNDER  INPUT  PARAMETERS  DELTA  AND  PARH(A). 
1-2,  THE  NUMBER  OP  FUNCTION  EVALUATIONS 
REQUIRED  DURING  THE  H0RK(5)  ITERATIONS. 
1-3.  THE  ESTIMATED  NUMBER  OF  SIGNIFICANT 
DIGITS  IN  OUTPUT  VECTOR  X. 

I-)l.  THE  FINAL  VALUE  OF  THE  MARQUAROT 

SCALING  PARAMETER  DESCRIBED  UNDER  PARM(l). 

ZEROS  AND  EXTREMA 
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1*5.  THE  NUMBER  OF  ITERATIONS  (I.E..  CHANCES 
TO  THE  X  VECTOR)  PERFORMED. 

SEE  PROGRAMMING  NOTES  FOR  DESCRIPTION  OF 
THE  LATTER  ELEMENTS  OF  WORK. 

INFER  -  AN  INTEGER  THAT  IS  SET,  ON  OUTPUT,  TO 

INDICATE  WHICH  CONVERGENCE  CRITERION  WAS 
SATISFIED. 

INFER  -  0  INDICATES  THAT  CONVERGENCE  FAILED. 
lER  GIVES  FURTHER  EXPLANATION. 

INFER  •  t  INDICATES  THAT  THE  FIRST  CRITERION 
WAS  SATISFIED. 

INFER  •  2  INDICATES  THAT  THE  SECOND  CRITERION 
WAS  SATISFIED. 

INFER  -  a  INDICATES  THAT  THE  THIRD  CRITERION 
WAS  SATISFIED. 

IF  MORE  THAN  ONE  OF  T»E  CONVERGENCE  CRITERIA 
WERE  SATISFIED  ON  THE  FINAL  ITERATION, 

INFER  CONTAINS  THE  CORRESPONDING  SUM. 

(E.C.,  INFER  •  3  IMPLIES  FIRST  AND  SECOND 
CRITERIA  SATISFIED  SIMULTANEOUSLY). 

ICR  >  ERROR  PARAMETER  (OUTPUT) 

TERMINAL  ERROR 

ieR»l29  IMPLIES  A  SINGULARITY  WAS  DETECTED 
IN  THE  JACOBIAN  AND  RECOVERY  FAILED. 
ICR«t30  IMPLIES  AT  LEAST  ONE  OF  M,  N,  lOPT, 
PARM(I),  OR  PARM(2)  WAS  SPECIFSCD 
INCORRECTLY. 

ICR- 132  IMPLIES  THAT  AFTER  A  SUCCESSFUL 
RECOVERY  FROM  A  SINGULAR  JACOBIAN.  THE 
VECTOR  X  HAS  CYCLED  BACK  TO  THE 
FIRST  SINGULARITY. 

ICR-133  IMPLIES  THAT  MAXFN  WAS  EXCEEDED. 

WARNING  ERROR 

IER-36  IMPLIES  THAT  THE  JACOBIAN  IS  ZERO. 

THE  SOLUTION  X  IS  A  STATIONARY  POINT. 
ICR-39  IMPLIES  THAT  THE  MARQUARDT 
PARAMETER  CXCCCDCO  PARH(3).  THIS 
USUALLY  MEANS  THAT  THE  RCQUCSTED 
ACCURACY  WAS  NOT  ACHIEVED. 

PRECISION  >  SINCLC  AND  DOUBLE 
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ZXSSQ  Is  a  finite  difference,  Levenber g-Marquardt  routine  for  solving 
nonlinear  least  squares  problems.  The  problem  Is  stated  as  follows: 


given  H  nonlinear  functions  f^ 

minimize 
over  X  ~ 

where  i  • 
be  estimated. 


f2,...,f^  of  a  vector  parameter  x. 

Is'a  vector  of  N  parameters  to 


When  fitting  a  nonlinear  model  to  data,  the  functions  f.  should 
be  defined  ae  follows:  * 

fl  it)  ■  1-1, 2,. ..,M  (l.e.,  the  realtjuals) 

where  Is  the  l^th  observation  of  the  dependent  variable. 

v^  •  tvJ.Vg . ''kv^  **  *  vector  containing  the  1-th  obser¬ 

vation  of  the  NV  independent  variables,  and 

g  la  the  function  defining  the  nonlinear  model. 

ZXSSQ  Is  based  on  a  modification  of  the  Levenberg-Harquardt  algorithm 
which  eliminates  the  need  for  explicit  derivatives. 

Let  x^  be  an  initial  estimate  of  x.  k  sequence  of  approximations 

m  a» 

to  the  minimum  point  is  generated  by 

.  •  nnnn  n«. 

where  J.  la  the  numerical  Jacobian  matrix  evaluated  at  x*^ 
n  • 

0^  la  a  diagonal  matrix  equal  to  the  diagonal  of 

for  lortao. 

la  a  positive  scaling  constant  (Narquardt  parameter) 

Ifhen  forward  differences  are  used,  the  Jacobian  Is  calculated  by 

J  1/2 

where  Uj  la  the  J-th  unit  vector  and  hj -max( | x j  |  ,0. l ) *eps 

(sec  scaling  comments  in  Programming  Note  2),  where  eps  la  the  relative 
precision  of  floating  point  arithmetic.  For  central  differences. 

la  used.  To  minimize  the  number  of  function  evaluations  required 
for  nonzero  lOPt,  a  rank  one  update  to  the  Jacobian  matrix  is  used 
where  approprlatet 

ZEIOS  SNO  EXTIENA 
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r' 


V,  •  •’n  *  TO 

where  <  -  x*'*'  -  x”. 

See  references: 

1.  Brown,  K.  H.,  and  Dennis.  E*.  "Derivative  free  analogues 

of  the  Levenberg'Har quardt  and  Gauss  algorithms  for  nonlinear 

least  squares  approximations" ,  Numerische  Mathematlk,  i8,  1972, 
289-297. 

2.  Brown,  K.  N.,  "Computer  oriented  methods  for  fitting  tabular 
data  in  the  linear  and  nonlinear  least  squares  sense".  Department 
of  Computer,  Information,  and  Control  Sciences,  TR  No.  72-13, 
University  of  Minnesota. 

3.  Levenberg,  K.,  "A  method  for  the  solution  of  certain  non-linear 

problems  in  least  squares".  Quart.  Appl.  Hath . ,  2,  I6il-l68. 

Harquardt,  D.  U.,  "An  algorithm  for  least-squares  estimation 
of  nonlinear  parameters",  S lAM,  1  1  (2)1963. 

Programming  Notes 

1.  Accuracy  in  evaluating  the  functions  f^  Is  critical  when  highly 

accurate  parameter  values,  Xy  are  required  (l.e.,  when  NSIC 

Is  greater  than  3  for  single  precision).  In  these  oases,  it 
la  advisable  to  evaluate  the  functions  in  precision  greater 
than  the  working  precision  (e.g.,  for  single  precision  ZXSSQ, 
double  precision  can  be  used). 

2.  In  reference  to  significant  digit  tests  for  the  x  vector  and 

the  sum  of  squares,  leading  zeroes  to  the  right  of  the  decimal 
point  are  counted.  For  example,  both  123.^5  and  0.001  23  have 
five  significant  digits.  Scaling  the  x  vector  in  the  functions 

f^  may  be  required  if  several  of  the  parameters  x^  are  much 

less  than  1.0  to  obtain  the  same  number  of  significant  digits 
for  each  of  the  Xj . 

3.  For  lOPT  4  0,  strict  descent  Is  enforced  by  increasing  the 

Harquardt  parameter  (see  PARM(I))  by  the  factor  PARM(2)^ 

until  a  smaller  sue  of  squares  la  obtained.  If  a  smaller  S3Q 
Is  obtained  immediately  during  Iteration,  la  decreased  by 

the  factor  PARM(2).  also  modified  by  the  ratio  of 

the  norm  of  the  gradient  2||J  f||  for  two  successive  iterations 
In  the  range  (PARH(2)'‘\PARH(2) ) ,  decreasing  as  the  gradient 

decreases,  and  increasing  ^3  the  gradient  increases. 
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.  For  IOPT-0,  a  formula  la  used  to  calculate  o  ,  and  no  check 

la  made  to  see  that  SSQ  has  been  reduced.  The  successive  Iterates, 

are  accepted  even  if  SSQ  increases.  IOPT-0  seems  to  be 

moat  effective  when  SSQ  has  a  minimum  near  zero.  In  other 
cases  and  when  IOPT-0  fails.  It  is  advisable  to  use  lOPT-i 
or  2. 


5.  The  Jacoljlan  has  typical  element  3f  /Jx.,  for  1-1, 2, ...,M  and 

*  j  2 

J-1,2,...,N.  The  Hessian  has  typical  element  3  ( SSQ  )  /  3Xj  3Xj^ , 
for  j  ,k-l  ,2  , .  . .  N.  liote  that  SSQ-f^f.  For  small  residuals 
f^,  the  output  matrix  XJTJ  Is  a  good  approximation  to  the  Hessian 
times  one  half. 


6.  For  normal  completion  (IER-0).  a  Jacobian  calculated  by  forward 
or  central  differences  (not  a  rank  one  update  version)  Is  returned 
In  ZJAC.  Also,  no  convergence  teats  are  performed  while  the 
rank  one  approximation  la  being  used. 


7.  If  a  terminal  error  other  than  IER-130  is  returned,  all  output 

parameters  contain  the  values  Indicated  for  the  Iteration  during 

which  the  terminal  error  occurred.  In  particular,  the  five 

values  returned  In  vector  WORK  should  be  inspected  whether 

a  terminal  error  was  returned  or  not.  If  lER«13i  la  returned, 

W0RK(3)  will  overeatlmate  the  number  of  significant  digits 

In  X  by  2  or  3  digits. 

«» 


a. 

9. 


Automatic  recovery  from  singularities  In  the  Jacobian  (l.e., 
an  entire  column  is  zero)  is  provided. 


Theoretically,  the  norm  of  the  gradient.  2||J  r)|,  should  be 
zero  at  the  least  squares  solution,  so  that  W0HK(1)  should 
be  inspected  before  a  solution  la  accepted.  Uaually,  values 


-it 

less  than  10  are  acceptable,  but  this  number  is  also  dependent 
on  the  scale  of  the  functions.  The  ersmabt  times  ona  half 
is  located  in  output  vector  U0RK(J*K}.  where  KoKAX((li«l  )*N/2.S) 


and  J-1 ,2, .. .,N.  WORK 


i(SSQ)/ix 


J* 


Example 


Suppose  the  nonlinear  model  g(Xiv^)  ••  X|*v^  «  sinlx^  vM  is  to  be 

fitted  to  5  data  points  y^  •  v^  ♦  sin(v‘)  •  (0.1)*(-1)^  where  v*-l. 
for  1-1 .2, . . . ,5.  Then  f  •  y,“f{x,v*)  and  for 

•  4  * 


ZXSSQ-6 
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APPENDIX  F.  DATA  PROCESSING  PROGRAMS 

The  programs  in  this  appendix,  TRANS. FOR  and  PRO.M,  are 
used  to  process  the  experimental  data  for  use  in  the  modeling 
program.  TRANS. FOR  is  a  FORTRAN  code  program  which  averages 
the  values  of  the  experimental  data  and  PRO.M  is  a  MATLAB 
program  which  normalizes  the  resultant  data,  from  TRANS. FOR, 
Co  ambient  and  converts  the  temperatures  from  degree 
Fahrenheit  to  degrees  Kelvin. 
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Program  TRANS 


c 

c  This  program  translates  data  into  free  format  and  smooths 

c  the  data  with  a  moving  window, 

c 

real*8  x(82) ,y(82) ,xl(61) ,yl(61) 
c 

open ( 8 , name= • pre3 . dat ' , status= ' old ' ) 
open ( 9 , name= • tvl . mat ' , status®  *  new • ) 
open ( 7 , name® • tv2 . mat ' ,  status® • new • ) 
open ( 1 1 , name® ' t v3 . ma t ' , status®  *  new ' ) 
open ( 10 , name® ' datam . dat • , status® ' new • ) 


c  Translate  the  Tvl  data 

do  i®l,61 
read(8, *)pl 
xl(i)=pl 
1020  enddo 

yl(l)®xl(l) 
yl(2)®xl(2) 
do  i=3,59 

yla)=0.2*(xl(J-2)+xl(i-l)^-xl(i)+xl(i+l)+xl{i+2)) 

enddo 

yl(60)»xl(6C; 

yj(6l)®xl(61) 

do  i-1,61 

write(9, I000)yl (i) 
c  write (10, ’k)  yi(i) 

enddo 

1000  format(F10.2) 

c  Translate  the  Tv2  data 

do  i-1,61 
read (8,*) pi 
xl(i)“pl 
1030  enddo 

yl(l)-xl(l) 
yi(2)«xl(2) 
do  i®3,59 

yl  ( i ) -0 . 2*  ( xl  ( i-2 ) +X1  ( i-l )  •►xl  ( i) +X1  ( i+i ) +xl :  i-t-2 ) ) 

enddo 

yl(60)®xl(60) 

yl(6l)«xl(6l) 

do  i®l,6l 

write(7,l0\0)yl(i) 
write(l0,*)  yl(i) 

enddo 
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for(Dat(Fl0.2) 
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c  Translate  the  Tv3  data 


do  i=l,61 
read (8,*) pi 
xl(i)=pl 
1040  enddo 

yl(l)»xl(l) 
yl(2)=xl(2) 
do  i=3 , 59 

y 1 ( i ) =0 . 2* ( xl ( i-2 )  +xl ( i-1) +xl ( i ) +xl ( i+1 ) +xl ( i+2 ) ) 

enddo 

yl(60)=xl(60) 

yl(61)=xl(61) 

do  i»l,61 

write (11, 1130) yl(i) 
write(10,*)  yl(i) 

enddo 

1130  fonnat(F10.2) 


close (8) 
close (9) 
close(lO) 
close (11) 

end 
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%  This  program  PRO.M  normalizes  the  smoothed  data  in  TV. MAT  and  TM.MAT 
%  to  the  ambient  value  which  is  the  first  temp  in  the  series  of  each 
%  and  then  converts  the  data  from  degrees  F  to  degrees  K. 

%  NOTE  -  THIS  PROGRAM  IS  RUN  AFTER  THE  TRANS. FOR  PROGRAM  IS  RUN 

load  tvl 
load  tv2 
load  tv3 
% 

tvsl=tvl-tvl  ( 1) ; 
tvsl=tvsl*0 . 535 ; 

% 

tvs2=tv2-tv2 (1) ; 
tvs2=tvs2*0 . 555 ; 

% 

tvs3=tv3-tv3 (1) ; 
tvs3=tvs3*0 . 555 ; 

% 

datam= [ tvsl *  tvs2 *  tvs3 ; 
datam 

save  datam.dat  datam  /ascii 
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